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I. INTRODUCTION 


Dynamic stall refers to the delay of stall onset beyond static stall angles 
experienced by airfoils and wings in unsteady motion. The phenomenon first drew 
serious attention in connection with helicopter aerodynamics when conventional 
methods of analysis proved incapable of accurately predicting performance for vehicles 
in high speed forward flight. The observed increase in overall lift could be explained if 
the lift on the blade moving opposite to flight direction was greater than predicted by 
steadv flow calculations. [Ref. 1] It was experimentally observed that the lifting 
characteristics of rotors could be adequately modeled by an airfoil oscillating in pitch. 
The basic mechanism is the shedding of a strong leading-edge vortex which distorts the 
pressure distribution as it travels over the upper surface of the airfoil and leads to the 
abrupt changes in lift and pitching moment usually associated with dynamic stall. 
[Rei 2] 

In addition to helicopter rotors, the dynamic stall phenomenon is observed in 
such diverse aerodynamic applications as wind turbines, jet engines, and rapidly 
maneuvering aircraft. Interest in expanding the design envelopes of fighter aircraft 
through both post stall maneuvers and extended conventional maneuverability has 
increased interest in dynamic stall research. Progress in this area depends upon 
improved knowledge of details of viscous flow over airfoils in dynamic stall. [Ref. 1] 
The basic airfoil dvnamic stall process is depicted in Figure 1.1, taken from Reference 
1. The flow-field pattern at anv instant is critically dependent upon the entire time 
history, so understanding the chronology of events including time prior to the vortex 
shedding is crucial. There are four basic stages: 

l. trailing-edge flow reversal, progressing forward along the airfoil, 


2. formation, growth, and shedding of the vortex, 


3. full stall, 
4. boundary layer reattachment. 
The specific characteristics of a given oscillating airfoil depend primarily on 


Mach number, 


= 


Revnolds number, 


G2 


reduced frequency, k, the ratio of the vertical velocity of the leading edge to the 
freestream velocity, k = oc/2U, where ® = circúlar frequency, c = airfoil 
chord, and U = freestream velocity, 
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Figure 1.1 Dynamic Stall Chronology (from Carr, Ref. 1). 


4. mean angle of attack, Uy, 


5. oscillation amplitude, 4,. [Refs. 1,3] 


While much of previous research has involved experiments, primarily on airfoils 
oscillating in pitch, advances in computer technology have increased the potential of 
computational methods to model details of unsteady flow fields accurately. Early 
analytic approaches based on inviscid models or semi-empirical analysis were very 
limited due to the complexity of the dynamic stall phenomenon and its interrelation 
with viscous effects. Moreover, boundary laver corrections are not appropriate due to 
the presence of large regions of separation. [Ref. 4] The Navier-Stokes equations can 
deal with flow separation and shock-boundary layer interaction, and the thin-layer 
approximation has been successfully applied to airfoils near maximum lift, including-- 
though somewhat less successfullv--conditions beyond maximum lift coefficient [Ref. 5]. 
The thin-laver Navier-Stokes equations may not be appropriate in the case of highly 
separated flows, however, and in the case of dynamic stall, the viscous layer is of the 
same order as the airfoil chord during the shedding of the strong leading-edge vortex. 
The full Navier-Stokes equations must therefore be solved. [Ref. 2] 

A Navier-Stokes problem solver has been developed by L. N. Sankar and his 
associates at the Georgia Institute of Technology and was made available for the 
present study. An implicit finite-difference procedure is used to solve the two- 
dimensional, Revnolds-averaged, compressible Navier-Stokes equations in strong 
conservation form. The flow solver features a moving, body-fitted coordinate system, 
an algebraic turbulence model (Baldwin-Lomax), and an alternating-direction-implicit 
(ADI) time-marching algorithm. It can also be used in an inviscid mode to solve the 
Euler equations. [Refs. 6,7] The goals of the present study were: 

1. Develop procedures for exercising the code on the Cray X/MP-48 computer to 


obtain output in a form suitable E E the pressure distribution with 
experimental results and for investigating details of the computed flow field. 


2. Po the code to steadv-state cases and to the conditions of existing dvnamic 
stall “experimental data [Ref. 8] and investigate the effect of varving input 
parameters. 

3; 


Obtain flow-field solutions under НЬ conditions of a series of wind tunnel 
experiments to be conducted at the, Fluid Mechanics Laboratorv of the 
National Aeronautics and Space Administration, Ames Research Center. 
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Il. MATHEMATICAL FORMULATION 


A. THE SOLUTION DOMAIN 

The first step in any numerical solution procedure is to discretize the physical 
domain. This involves the selection of suitable boundaries and the choice of 
appropriate nodal points within the boundaries to describe the characteristics of the 
physical domain accurately. In the solver used here a solution grid is constructed 
about the airfoil, which is then rotated along with the airfoil for dynamic calculations. 
With the solution domain defined, it is then necessary to specify conditions at the 
boundaries. For time dependent problems such as those considered here, initial 
conditions at the nodal points must also be defined. 

1. Algebraic Grid Generation 

A mathematical transformation from the physical plane to a computational 
plane introduces significant advantages for numerical solution. Chief among these is 
the fact that the physical boundaries are mapped into rectangular surfaces in the 
transformed plane. Grid points can be concentrated in physical regions experiencing 
the largest gradients while maintaining the simplicity of uniform spacing in the 
computational plane. [Refs. 9,10: pp. 519-520] Proper design of the grid 1s crucial to a 
stable, accurate solution. The version of Sankar's program used for this study emplovs 
an algebraically-generated C-grid. This grid is quickly produced, easily rotated by 
simple sine/cosine relationships with the angle of attack, and allows a high degree of 
clustering in the normal direction to cover adequately the thin boundary laver of high 
Reynolds number flows. 

The procedure used generates a sheared parabolic coordinate system based on 
an input airfoil shape. Two points are first found on the airfoil: point N halfway 
between the nose and the center of curvature of the nose, and point T at the trailing 
edge. The cut along the airfoil wake is chosen to be tangent to the mean camber line 
at the trailing edge. The airfoil shape and wake are then unwrapped to form a surface 


S(X) in an intermediate X-Y plane using the following transformation: 


ШЕ O) = AZ ZN (eqn 2.1) 


ІШ 


where z = x + 1y locates a point in the physical (x-v) plane. This surface will then lie 
along the X-axis with a shallow bump at the origin which coincides with the singular 
point zx. Figure 2.1 illustrates this mapping for a symmetric airfoil. Cubic 
interpolation is used to identify additional points and smooth the surface S(X). Lines 
parallel to the Y-axis and lines equidistant from the surface S(X) are then constructed 
before the sheared cartesian grid is mapped back to the physical plane. Each point in 
the physical domain is assigned to a corresponding point in the computational (€-n) 
plane through the intermediate plane relationship. The intermediate and 
computational planes are illustrated in Figure 2.2. Unit spacing is assumed in the 
computational plane, simplifving calculations. For solutions of viscous flows, the § = 
constant lines are retained while the n = constant lines in the physical plane are 
redistributed to place points within the boundary layer. The first point off the solid 
surface is placed at a distance specified by the user and the remaining lines are 
exponentially distributed to the far-field boundary. [Refs. 6,7: pp. 40-44] 

The final mesh used is shown in Figure 2.3. The number of points defined in 
the 3-direction is 161 of which 159 are used for calculations, and 41 points are defined 
in the y-direction. A detail of the grid around the airfoil is shown in Figure 2.4. The 
spacing along Constant lines is relatively fine near the leading edge, where the 
sharpest gradients are expected, and coarse at the trailing edge. In order to perform 
calculations in the computational plane, a transformation 15 required to map each 
point to the corresponding point in the physical plane. In the present solver, a 
numerical approach is taken: finite differences are used to relate the mesh spacings in 
the two planes through transformation metrics. This method 1s outlined in later 
sections. 

2. Initial Conditions and Boundary Conditions 

The initial conditions for steady-state calculations are assumed to be the 
freestream conditions, and inaccuracies introduced by this crude starting approximation 
are removed by viscous effects after the airfoil is impulsively started. (Inviscid 
calculations require proper boundary conditions and artificial dissipation to correct the 
solution). [Ref. 7: p. 13] In the case of dynamic calculations, the initial conditions are 
obtained from an asymptotically converged steady-state solution at the minimum 
airfoil angle of incidence. The solution is then marched in time as the oscillation cycle 
begins from this angle. Boundary conditions are explicitly applied at each time step on 


the solid boundary, at the far-field boundary, and in the airfoil wake, where the grid 
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Figure 2.1 Symmetric Airfoil in the Phvsical Plane (top) 


Unwrapped in the Intermediate Plane (bottom). 
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Figure 2.2 Construction of the Grid in the Intermediate Plane (top) 


Corresponding Computational Plane Grid (bottom). 
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Figure 2.3 Al 
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generation procedure introduced a “cut” in the physical domain between the airfoil 


trailing edge and the downstream boundary. (See Fig. 2.5) [Refs. 6,11] 


B. GOVERNING EQUATIONS 
The unsteadv, compressible, two-dimensional Navier-Stokes equations are written 


in cartesian coordinates as: 


Ó 
ES. pp 
ct 


P = | pu ir pv |, 


2 
pu“ + pu puv-t,, 


+ p-T 
цоог ШОШЫП xy ЧЕ Q. (Est рит, E Q, 


puv-t, 


2 
| 
Зэ 


е (eqn 2.2) 


where p, u, v, and e are density, velocity components, and total energv per unit 


volume; T and со, are components of the shear-stress tensor derived by 


t 
хх’ уу y 
application of Stokes’ hypothesis; and Q, and Q, are the heat-flux vector components 


given by: 





Ч! ES 
EU (eqn 2.3) 


in which M is the coefficient of viscosity, Pr is the Prandtl number, and a is the speed of 
sound. [Refs. 12,10: pp. 480-481] The specific heat ratio, y, is assumed to be 1.4. 


Expanding and rearranging Equation 2.2 vields [Refs. 6,7: pp. 9-10]: 


д4 вое 2,6 = 0 R ІН 0,5 (eqn 2.4) 
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with 


а = [р Е = [ри С = [ру | 
pu pu`+p puv | 
ру puv ptp | 
е (ESP е+р)у| 

В hi 1 > = [о 

| Ux Us 
Kä т, 
К, 5, = 


and 
O Vp + Кд (а) 
о Kô (a?) 
1 
К = iere 
(y-1)Pr (eqn 2 


In this form F and R contain the flux and viscous stress terms in the x-direction and G 
and S the flux and viscous stress terms in the v-direction. Note that the viscous terms 
are contained on the right-hand side (setting them equal to zero produces the Euler 
equations). The governing equations are given in strong conservation form. That is, 
the coefficients of the derivative terms are constant, and the equations may be written 
in a form expressing the divergence of the physical quantities of mass, momentum, and 
energv (Eqn. 2.2). Such equations may be expressed in a corresponding integral form 


through the divergence theorem: 
J,V * Pdv 2 0 = J.P ends (eqn 2.6) 


The differential statement of the governing equations is not valid across discontinuities 
(1e., shocks); however the integral expression is. [Ref. 13: pp. 406-408] Using the 
conservative form means that the finite-difference approximation will ensure 
conservation of the physical properties and thus satisfy the corresponding integral form 
of the governing equations. [Refs. 14,10: pp. 50-52] This gives weak solutions (solutions 
of partial differential equations that include discontinuities) [Ref. 10: pp. 139-141] in 
the vicinity of shocks. Such methods are called shock capturing, and their use 15 


important in treatment of the compressible dynanuc stall problem. 
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AS discussed earlier, considerable simplification of the numerical solution mav be 
realized by employing a transformation into a rectangular plane. It can be shown that 
the strong conservation-law form may be restored following such a transformauon 


[Refs. 14,9: p. 254]. The general transformation is given by: [Ref. 6] 


5 = S(x.y,t) 
Ц = ЧУ 
c M (eqn 2.7) 


and the Jacobian and metrics of the transformation by: 


J = Ny Пу = l/eyq — xqvg) (eqn 2.8) 
Sx = I¥y 

n 

5, E Хх т Ve, 

ШЕЕ УЕ 

We Jee 

get, yn, (eqn 2.9) 


Here the subscripts indicate partial differentiation. with respect to the subscripted 
variable. The Jacobian represents the magnification factor of mesh elements between 
the physical and computational planes [Ref. 10: p. 530]. The time- and spatial- 


derivatives of any quantity q are given by: 


Ч, = 9; + 95, + qq 
9; = 95, + а, (сап 2.10) 


Ihe governing equation (Eqn. 2.4) in the transformed plane may then be expressed as: 


x Ж ж - Ж Жж 
9.4“ + ôg F + On G CR + Cu (еап 2.11) 


where 

17-4: 

FY = (& F + EG Say 
6* = (И.Р + 19,6 +10) 


R S 
... 


SUM 


(Š R + & S)/J 


| 


ж 3 
S (HR + n,5)/J (eqn 2.12) 
This constitutes a coupled. highly non-linear system of equations. The flux vectors F* 
and G* may be stated in a less complex form bv introducing the contravariant velocity, 
which is required to compensate for grid motion. The contravariant velocity 


components are given bv [Ref. 9]: 


| 
IW 


С 4 do S. u ps SV x S (u—x.) + S (T) 
MEANS AU Is V Um SE UL (VE) (сап 215! 


Using these velocities, F* and G* тау be rewritten, as used in їе р m 


implementation: 
1 de l D 
ES ча puU + GP С* = 3 puV * p 
РУС + & p PVV t "n.p 
Cae Ир (еруу кш (eqn 2.14) 


C. NUMERICAL SOLUTION 

The finite-difference scheme emploved in Sankar’s problem solver is based on the 
Beam-Warming approximate factorization algorithm ([Ref. 10: pp. 489-496] as 
implemented by Steger [Ref. 9]. A theoretical description is given in References 6 and 
7. Solving the governing equation in the transformed plane (Eqn. 2.11) requires 
determination of the metrics and Jacobian, relating the variables to corresponding 
quantities in the physical plane through Equations 2.12 and 2.14. This is accomplished 


through Equations 2.8 and 2.9, where central difference approximations are used to 


compute the spatial derivatives ХЕ, УЕ, Хи, апа yq at each interior grid point, and x, 


and v, are the grid velocity components. One-sided differences are used at the 


Г 
boundaries. The governing equation (Eqn. 2.11) 1s parabolic (hvperbolic when the 
viscous terms are neglected) [Ref. 10: p. 139], lending itself naturally to an implicit 
solution approach [Ref. 14]. This allows the use of a larger time step than explicit 
schemes, and the requirements on step size are generallv set bv accuracy rather than 
stability requirements [Refs. 2,9]. The use of explicit boundary conditions in the 
present implementation introduces a step-size stabilitv requirement. Given an assumed 
flow field at time tp an alternating-direction-implicit (ADI) procedure is used to 
advance the solution to time level t_ , ¡. Elements of the vector q* are then given in 
terms of values at time t, and increments, ie. газа 1 qu + (Agent). In 
Sankar's program, the approximation method used is an Euler implicit finite-difference 
scheme [Ref. 10: p. 98], with central spatial differences and backward time differences 
for the derivatives in the governing equation. The method is thus first-order time 
accurate and second-order accurate in space. [Ref. 7: pp. 21-23] 
1. Linearization and Factorization 

Inspection of the terms of the flux and viscous stress vectors of the governing 
equation shows them to be very non-linear functions of the unknown vector q*. The 
solution procedure involves the use of truncated Taylor series expansions to linearize 
the equation [Ref. 10: pp. 490-491]. The viscous terms, however, are treated explicitly. 
evaluated from the solution at the previous time level, rather than following Steger’s 
fully implicit treatment [Ref. 9]. This leads to improved computational efficiency but 
requires implicit smoothing for stability at moderate to high Reynolds numbers 
[Ref. 6]. The Euler implicit form may be obtained from Taylor series expansion at time 
t; p Replacing the partial time derivative with the backward difference operator V, 


vields: 
(Aq*) ^! 2 (qs? ^L — (q^ - AqV. (q^ * ! + O(At)? (eqn 2.15) 


Substituting Equation 2.11 into Equation 2.15, after replacing the spatial derivatives 
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with the central difference operators бе апа бұ and the time derivative with the 
backward difference operator. produces the governing equation in non-linear, time- 
differenced form [Refs. 7,10: pp. 22-23,490-491]: 


44911 = — Aw 6c fpein t! got (Саул!) 
- Ad (ROT! + eS + O(At) (eqn 2.16) 


Tavlor series expansions of the flux vectors are now introduced: 


(келті = (F°)" + ела + (АО 
G = (С) + o C C O(AD (eqn 2.17) 


Here [e EP and AGT are Jacobian matrices, evaluated under the assumption of a 
perfect gas [Refs. 7,10: pp. $7-88,491], which will be denoted [A] and [B] respectively. 
Substituting Equation 2.17 into 2.16 and placing the implicit terms on the left side and 


explicit terms on the right side gives the linearized system 


(1+ АА] + 46 В} (Ад = В (eqn 2.18) 


where 
RP = — Аҳд (Е) + o da Add (RENT! 1 пе) 


Since, as mentioned previously, the viscous terms are treated explicitly, they have been 
grouped together in {R}" with the other terms to be evaluated at the known time level 
г. [№е[5. 6.7: РВ. 23-23] 

Although linearized, the svstem of equations is now in block pentadiagonal 
form and still computationally expensive to solve. Approximate factorization is used to 
write the implicit scheme as a product of one-dimensional factors so that the problem 
may be solved by a sequence of relatively simple operations, subject to the additional 
error introduced by the factorization approximation. [Ref. 14] Factoring the left side 


of Equation 2.8 then gives [Rel | 
({I] + Atée{A]) (0 + Até,[B]) (Aq*}" 7! = (RJ (eqn 2.19) 
The Beam-Warming algorithm is now implemented in three steps [Ref. 10: p. 494]: 


25 


GEO 


(Ш + АА} (49**} = (К) (eqn 2.20) 
Sep 2 

(Ш + Ах БИА”) АЫ (да+) (eqn 2.21) 
Step 3 

Я а" (eqn 2.22 


Notice that steps one and two represent sweeps in alternating ¢- and n- directions, 
leading to the name alternating-direction-umplicit for this method. Since the difference 
operators б апда бү represent central difference approximations, they lead to svstems 
of block tridiagonal matrix equations composed of 4X 4 submatrices, which are easily 
solved. 
2. Application of Boundary Conditions 

On the solid surface the flow tangencv condition applies for the Euler 
equations, and in addition, the no-slip condition applies for the viscous flows 
considered here. This means the fluid and solid surface must have the same velocity at 
the common boundary, which is to say the contravariant velocity components U and V 
(Eqn. 2.13) are both zero. The surface is assumed adiabatic so Q e = 0. It is also 
assumed, as reasonable approximations for high Reynolds number flows, that dp -0 
and др = Q at the surface. Pressure and density at the surface are then numerically 
determined by two-point extrapolation, E (4p; — Р; 32/3 апа р: = (3p; » - 
P; 3)/3 internal energy is calculated using the relation p ^ (y — 1X (e — 0.5p(u* + ү2)). 
Boundary conditions are specified following each ADI sweep for the incremental 
quantities {Aq*} and (Aq**) by setting them equal to zero on the solid surface. 
[Refs. 6,7,11] The far-field boundary conditions are undisturbed freestream conditions. 
Conditions along the cut in the C-grid are found by averaging the solution at the 
Nearest interior points on either side of the cut. [Ref. 7: pp. 27-29] The effect of this 
scheme is to update conditions at the boundaries explicitly. The explicit boundary 
conditions are used for their simplicity in the code, although implicit conditions would 


generally be preferable for stability and accuracy of the numerical solution [Refs. 9,15]. 
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3. Artificial Dissipation 

When Tavlor series expansions are substituted into the finite-difference 
approximation of a partial differential equation, and the terms are rearranged to 
produce the original differential equation plus a truncation error, the lowest order 
truncation error term may contain a second derivative that makes the term similar in 
appearance to the viscous terms in flow equations. Thus the use of finite-difference 
approximations introduces an “artificial viscosity” into the solution. [Ref. 10: pp. 
89-92] One-sided difference schemes (upwind differencing) are commonly used to 
produce this implicit dissipation. Additionally, terms may be added explicitly to 
suppress oscillatory behavior of the numerical solution. Central difference 
approximations, given by (u.4, — u._,)/2Ax, tend to decouple the odd from the even 
terms, since only every other term is used in the difference formula. For example, the 
sequence of nodal values [—1,+1,—1,+1,..+ would give а central ddl NE 
approximation to the first derivative of zero. [Ref. 14] Solutions of the Navier-Stokes 
equations may “blow up” due to numerical oscillations if the mesh is not fine enough in 
areas of large pressure gradients. lo suppress high frequency numerical oscillations, 
fourth-order explicit dissipation terms are commonly added to algorithms. [Ref. 10: pp. 
105, 486] Physical dissipation diffuses energy, and artificial dissipation reduces gradients 
in the flow-field solution whether such diffusion is physically correct or numerically 
induced. The explicit method of adding dissipation has the advantage over one-sided 
differencing of making the physical approximations clearer and permitting some control 
over the amount of non-physical (numerical) dissipation. [Ref. 4] Although viscous 
terms provide a dissipative mechanism, some additional numerical dissipation is always 
necessary; the degree should be kept to the minimum required. [Refs. 14,10: p. 92} 

The Navier-Stokes code compiled by Sankar follows Steger in the use of 
artificial dissipation, with changes to prevent overshooting in the vicinity of shocks. 
[Ref. 16] Dissipation is employed for four main purposes, the first three of which have 
been mentioned previously: 


D М press high frequency oscillations in the solution caused by the use of central 
ifferences. 


2 CHE for the incorrect initial conditions, especially when solving inviscid 
OWS. 

3. Allow explicit treatment of viscous terms. 

4. Alleviate restrictive stability bounds on the use of explicit damping. [Ref. 9] 


„ә 
<> 


The first two reasons pertain to the right-hand (explicit) side of Equations 2.19 and 
2.20 and the second two reasons to the left-hand (implicit) side of Equations 2.19, 2.20, 
and 2.21. The explicit dissipation used is a combination of second- and fourth-order 
terms. Fourth-order smoothing is normally used for accuracy, since second-order 
terms tend to smear rapid pressure variations near the leading edge. In large pressure 
gradients such as shocks, however, these terms lead to overshooting in the pressure 
distribution. To prevent this, the second derivative of pressure 1s used to control the 
degree of dissipation, turning on the second-order term and suppressing the fourth- 
order term in the vicinity of shocks [Ref. 6]. 

Upwind difference schemes may be rewritten as the sum of central difference 


approximiations to the first and second derivatives as follows: 
(u.—u,_)/Ax = (u4,—u_))/2Ax—(u4y, + 2u,—u_,)/2Ax (боп 225) 


The second term on the right-hand side may then be regarded as an implicit dissipation 
term, giving a central difference scheme the artificial viscosity characteristics of upwind 
differencing. Second-order implicit smoothing terms are added to the left-hand side of 
Equations 2.19, 2.20, and 2.21, both to allow the use of explicitly evaluated viscous 
terms and to permit the use of larger magnitude explicit damping terms without loss of 
stability. The dissipation terms are of at most the same order as the truncation errors 
of the finite-difference approximation involved [Refs. 6,7: p. 30-33]. The coefficients of 
both implicit and explicit dissipation terms, £j and £p, are user selectable options that 


control the magnitude of the artificial viscositv. 


D. TURBULENCE MODEL 

Ihe Revnolds form of the Navier-Stokes equations introduces new unknowns in 
the form of turbulent stress and heat-flux terms, requiring additional relations to make 
solution of the system possible. This is known as the closure problem, and the 
approach usually taken to resolve it is rurbulence modeling. [Ref. 10: pp. 207-208, 
221-235] Under the Boussinesq assumption, the coefficients of viscosity and thermal 
conductivity are replaced by the combinations Н + H, and k + k, representing the sum 
of laminar and turbulent components. In practice the thermal conductivity is usually 
determined from the relationships k = с,Д/Рг and k, — cyIt/ Pr. In the solver used 


here a constant value of Pr = | is assumed, and total viscosity is used in Equation 2.3 
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rather than directlv calculating the conductivitv. For the present calculations, the 
viscous work and conductivity terms R, and 5, (Eqn. 2.5) were set equal to zero. The 
laminar viscosity is easily found from the Reynolds number and freestream velocity, 
and an algebraic turbulence model is used to determine the turbulent eddy viscosity. 
Turbulence modeling lies at the frontier of research in computational fluid 
dynamics. The use of Navier-Stokes solvers as practical, predictive tools depends on 
the development of adequate turbulence models. For many uses, however, the 
dissipative, diffusive mature of turbulence has enabled simple algebraic models to 
perform bevond their theoretical limitations [Ref. 17]. In the present case the complex, 
history-dependent viscous effects dominating the flow are clearly beyond the capacity 
of simple algebraic models based strictly on local flow properties. Unfortunately, there 
are no completely satisfactory alternatives. The Baldwin-Lomax model was chosen for 
simplicitv, while acknowledging that it might be unsuitable for the massively separated 
flows experienced in dynamic stall [Ref. 6]. When it was discovered that the maximum 
lift was underpredicted, a modification was incorporated which will be discussed 
following a description of the original Baldwin-Lomax model. 
1. Baldwin-Lomax Model 
The Baldwin-Lomax model is based on Cebeci’s two-layer model, with H 
given by the inner-laver formulation out to the smallest normal distance at which the 
inner- and outer-laver formulas give the same value for the viscosity. In the inner layer 
a simple mixing length formula is used with the eddy viscosity proportional to the 
magnitude of the local vorticity times the square of the normal distance from the 


surface: 


(ит): ег = РЕЈ 


where the mixing length 
{ = кур (едп 2.24) 


and D is the Van Driest damping function, given by: 


D = [1 — exp(-y * JA* J] (eqn 2.25) 


(22 
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The damping constant AT = 26, and the von Karman constant, К, 15 taken as 0.4. In 


the outer layer, the locally constant eddy viscosity 1s given by: 
= > 7 2315 
(MT)outer KC PF wakel Klep() (сап 2.26) 
The key feature is the definition of a function 
F 


: 
= 1 , - 10:25 > 
wake ^ 2™N(YmaxF max: Gwk¥max" ait /F max) (eqn 2.27) 


where ээ is the maximum of the function 


F(y) = yla|D (eqn 


tJ 


28) 


and Yio, is the value of y at which F,,,, occurs. The Klebanoff intermitancy 


ax 
function Fg. has the effect of causing the eddy viscosity to approach zero in the far 
field. Сағ is the difference between the maximum and minimum velocities in the 
velocity profile. The values of the constants used are K = 0.0168, Cop = 1.6, and 
ШО - 92 Refs. 2.6.7: pp. 16-19] 

Modified Model 


It has been observed that algebraic turbulence models of the Cebeci-Sniuth 


г 


type [Ref. 10: p. 225] perform well in attached flows at moderate angles of attack but 
overpredict both the rise in turbulent shear stress in adverse pressure gradients, and the 
stress decrease when pressure gradients are relieved [Ref. 18]. The Baldwin-Lomax 
model was found by Sankar to cause an early prediction of flow separation at high 
angles of attack. It was assumed that this was due to underprediction of the turbulent 
length scale at high angles. The function F, defined in Equation 2.28, has in some 
cases been found to possess a number of local maxima, leading to large length scale 
variations depending on the maximum chosen. This presented a possibility for simple 
modification to increase the outer-laver velocity and length scales. The definition of 
Fax Was redefined to be the value of the function F when the product yF(y) 1s 
maximum. In preliminary investigation this change was found to increase stall angles 
but to lead to premature reattachment during the downstroke in dynamic calculations. 


For this reason its use is recommended only until stall onset. [Ref. 19] 


ES 


Ill. DESCRIPTION OF THE CODE 


The version of Sankar’s code used in the present study is presented in Appendix 
A. Appendix B contains instructions and guidelines for running the program. 
Modifications were made in the output scheme to produce output of pressure 
coefficients and flow-field plotting data at a specified number of equal time intervals 
during each cycle, and job cards were prepared for save, restart, and various output 
options. Other features of the program were not altered, although comments were 
added. Subroutines in the program supplied by Sankar were fully vectorized for the 
Cray Computer with the exception of subroutine EDDY. Execution times =omeene 
NASA Ames Research Center's Cray X-MP/48 computer, for full viscous dvnamic 
calculations, are approximately 0.318 seconds per time step. 

All program variables are nondimensionalized. This has the effect of multiplving 
the right side of the governing Equation 2.11 by a factor of Re'! but does not affect 
the basic form of the equation [Ref. 10: pp. 191-193]. The nondimensionalization 15 
effected as follows: 

density with respect to P-o 
2 length mith respect toe 
3. time With respect to c/ag 
J. velocities with respect to ago 
3. total energy with respect to fo aco? 
Under this scheme the number of time steps required to complete an oscillation cvcle 15 
given бу л/(к х Мо * At). 


A. MAIN PROGRAM 

The program begins execution by reading the input data. This consists of the 
dynamic stall parameters, Ug, Uy. k, M, and Re; program control information, including 
time-step size, explicit viscosity coefficient, and flags for features such as the modified 
Baldwin-Lomax turbulence model; and the airfoil coordinates and related information. 
Table 1 summarizes the subroutine calls. The grid generation routine, AIRFOL, 1s 
called, and for viscous flows, CLUSTR is then called to recompute the constant-n grid 
lines so that the first point from the airfoil surface is at the user-specified distance. The 


initial conditions are established next, and if the program is being restarted from 
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previous calculations, the values of time and the unknown flow-field vector q* are 
replaced with those now read from logical unit 7. Finally, subroutine METRIC is 
called to compute the metrics and transformation Jacobian, and the iterative 


calculations begin. 


TABLE | 
SUBROUTINE CALLS FROM THE MAIN PROGRAM 


Program MAIN 
AIRFOL generate grid 
Cres tk cluster grid points near the surface for viscous calculations 


METRIC compute metrics and Jacobian 


Loop: 
ROTGRID rotate grid to current angle of attack 


METRIC recompute metrics 
SLPS perform ADI sweeps 
WALLBC applv boundary conditions on surface and cut 


ММ ОЛ output surface pressure distribution 


LOAD compute integrated lift, moment, and drag coefficients 





The iterative solution loop is designed to execute the ADI sweeps for a number 
of steps specified in the input data. For restarts of dynamic calculations, the iterative 
solution loop begins its first cycle by rotating the grid (subroutine ROTGRID) based 
on the time read from the previously saved solution and the frequency of oscillation. 
No rotation is required at the initiation of dynamic calculations or for steady-state 
calculations since the initial conditions set the freestream flow direction at an angle 
equal to the initial angle of attack. On each subsequent iteration the grid is rotated an 
incremental amount based on the size of the time step. Following each grid rotation, 
subroutine METRIC is again called to recompute the metrics since the computational 
plane is stationary. The actual solution is now accomplished by a call to SLPS, and 


boundary conditions are applied at each step by WALLBC. The completed solution 


may now be used to produce normal program output, at intervals specified in the 
program. CPPLOT is called to output the surface pressure distribution, and LOAD 
computes the integrated lift. drag, and moment coefficients. After exiting the loop, 
various output options may be selected, including the option of saving the current 
solution for a subsequent restart; creating the plotting files; and writing the complete 
flow-field solution, including the velocitv magnitude, eddy viscositv, and normal 
distance from the surface at each grid point. The present version of the program is 
modified to exit at a number of equal time intervals each cycle, as specified within the 
program. This was done to generate data files for PLOT3D (an Ames Research Center 
plotting code) at equal phase intervals, and also to save the current solution so it may 


be written to Cray tapes for possible future use. 


B. GRID GENERATION 

The call to subroutine AIRFOL initiates the basic grid generation process. 
AIRFOL begins by reading the airfoil shape input parameters. These include the 
number of points on the upper and lower surfaces and a flag indicating airfoil 
synimetry. For svmmetric airfoils, only the upper surface coordinates are read, and the 
lower surface coordinates are defined as (XL, YL) = (XU, —YU). Еа 
coordinates are scaled with respect to the airfoil chord length, c, and the actual grid 
generation begins. The sequence of subroutine calls is contained in Table 2. 

The first step in the grid generation is to fill in the definition of the airfoil 
surface. Subroutine SING is called to determine the singular point (xy, yy), defined 
earlier as lying midway between the nose and the origin of the leading-edge radius. 
SING applies a three-point parabolic curve fit to the leading-edge points to determine 
the radius of curvature. The angle of the trailing-edge upper surface relative to the 
leading edge and the average of upper-surface leading- and trailing-edge slopes are also 
computed for use in determining the wake shape. AIRFOL next calls subroutine 
TABINT to compute the additional points required. TABINT first. performs the 
mapping to an intermediate plane described by Equation 2.1. The distance required 
between points to place 97 points on the airfoil surface is calculated. Then the actual 
determination of the points is performed by a standard polynomial interpolation 
routine, TAINT. Finally, TABINT maps the smoothed surface back to the physical 


plane. 


TABLE 2 
GPID GESNPEFRATISO SSSDPROUTINESCALLED 


AIRFOL - called by MAIN 
SING - determine singular point (Xv, vx.) 


TABINT - find additional points on airfoil surface 
TAINT - interpolation routine 
WRAP - calculate stretching in N-direction 


РОТЕ called bv MAIN for viscous flows 
STRTCH calculate new stretching factor 


TAINT locate new normal grid point locations 





A smooth airfoil shape has now been defined. [t remains only to determine the 
shape of the wake--the “cut” in the physical domain--before constructing the grid. 
AIRFOL uses the trailing-edge angle computed by SING to calculate a wake shape 
allowing the flow to leave the trailing edge smoothly. This cut generally corresponds 
to the tangent of the mean camber line, and for symmetric airfoils, it is just the 
extended chord line. [Ref. 7: p. 41] It remains at a fixed angle, and the entire grid is 
rotated along with the airfoil. (See Fig. 2.4) 

Construction of the grid is finally completed by subroutine WRAP. The 
stretching in the normal direction is computed. Then the airfoil and wake are 
unwrapped using Equation 2.1. This creates a grid consisting of constant-n lines 
equidistant from the unwrapped surface and constant-§ lines parallel to the y-axis. 
This grid is used for inviscid (Euler) calculations. Following the return from WRAP, 
the last step performed by AIRFOL is to map the entire grid back into the x-y plane, 
and the coordinate axis is shifted to the quarter-chord point. 

If the Reynolds number is greater than zero, program MAIN calls subroutine 
CLUSTR to construct new constant-y lines. (Setting the Reynolds number to zero or 
a negative value is a flag for inviscid calculations). For each grid x-coordinate, the old 
Stretching in the normal direction is computed; that is, the physical distance from the 


surface to each point along constant-G lines is calculated. Then subroutine STRTCH 
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is called to calculate the new stretching. In STRTCH the user-input distance of the 
first point off the wall and the distance just calculated bv CLUSTR to the farthest grid 
point are used as the inner and outer bounds of the grid. A stretching factor R is 
desired to give a geometric progression of grid spacings such that, given the distance of 
the first point off the wall, the ratio between successive spacings is R. This correct 
stretching factor 1s calculated by iteration using Newton's method. Once the stretching 
factor has been determined, CLUSTR calls the interpolation routine TAINT to locate 
the new grid coordinates. 

Construction of the grid in both the physical and computational planes is now 
complete. The phvsical grid is easily rotated by subroutine ROTGRID using the 
relations x = x cos (Ad) — y sin (Adu) and y = v cos (Ag) + x sln (AG). The only 
step remaining is the calculation. of the transformation metrics and Jacobian. 
Subroutine METRIC is called by MAIN after construction of the grid and after each 
rotation. METRIC performs calculations as outlined in the previous chapter under 
“Numerical Solutions”. The grid velocity components are x, and y,, determined from 
the angular velocity and grid coordinates in the physical plane (since coordinates were 
defined with respect to the quarter chord). Central differences are used at interior 
points and one-sided differences at the boundaries--two-point downstream and three- 
point elsewhere--to compute ХЕ. Хүр УЕ, апа Уп: This calculation 1s simplified by the 
assumption of a unit grid spacing in the computational plane. Finally, the relations 


given bv Equations 2.8 and 2.9 are used to compute the Jacobian and metrics. 


С. COMPUTATIONAL STEPS 

The basic features of the computational scheme Were outlined in the preceding 
chapter. The organization of subroutines to implement this scheme is presented in 
Table 3. Program MAIN calls subroutine SLPS to perform the primary calculation 
steps. SLPS controls the execution of the ADI sweeps, assembles the matrices, collects 
the damping terms, and calls the matrix solvers. The value of the implicit damping 
coefficient is set in relation to the explicit coefficient. In the present implementation £j 
= 20=-. If the reduced frequency is less than .001, a local time-step option is used to 
accelerate convergence to a steady-state solution by basing the time step on the local 
flow conditions. The calculations begin by calling subroutine DISSIP to compute the 
explicit dissipation terms described in the last chapter and add them to the right side of 


the governing equation (Eqn. 2.20). In the -direction DISSIP uses a blend of second- 


and fourth-order terms. А switching function is computed, based on the second 
derivative of pressure. This function causes the fourth-order terms to be set to zero 
when the pressure gradient exceeds a value specified within the program. Otherwise 
the second-order terms are set to zero. The large pressure gradients experienced in the 
vicinity of shocks are not expected in the n-direction, so only fourth-order terms are 


added. The arravs DQ1-DQ4 contain the combined dissipation terms. 


TABLE 3 
ЭШЕ БӘ ЕЮ В COMPUTATIONAL STEPS 


SIE PS - called by MAIN 
DISSIP - compute explicit dissipation terms 
RESI - compute inviscid right-hand (explicit) terms 
STRESS - compute viscous right-hand terms 
EDDY - determine viscosity coefficient 


AMATI - compute Jacobian matrix [A] 


МАТЕХІ  - invert assembled matrix in the -direction 


AMAT?2 - compute Jacobian matrix [B] 


МАТКХ2 invert assembled matrix in the y-direction 


WALLBC called by MAIN to apply surface boundary conditions 





Subroutine SLPS next calls RESI to compute the inviscid right-hand terms at the 
known time level (Eqn. 2.18). Working first with the £-derivative flux terms, RESI 
begins by calculating the contravariant velocity component U, in the computational 
plane (Eqn. 2.13). The vector F* (Eqn. 2.14) is then formed using the contravariant 
velocity. Then АЦО: (Е*)"] is Computed using standard central differences and added 
to the arrays DQ1-DQ4. These steps are repeated in the n-direction, using the 
contravariant velocity component V, and Atle (G*)"] is added to the DQ arrays. 

For viscous flows, subroutine STRESS is called next to compute the viscous 
terms and add them to the right-hand side of the equation. The first step in STRESS 
is a call to EDDY to compute the total (turbulent and laminar) viscosity. (EDDY is 


not called on the first ten steps of the initial dynamic run to allow the solution to settle 
a bit). The eddy viscosity routine begins by initializing the viscosity throughout the 
field to the laminar value calculated from the Mach number and Reynolds number. ЇГ 
laminar calculations are desired, the subroutine should be exited at this point. A 
completely turbulent flow field is assumed, so no transition point must be determined. 
The subroutine works with one grid $-location at a time. The calculations performed 
were described in the last chapter under “Turbulence Model”. The variables defined in 
that section are initialized to low values which will be reset as the flow field is scanned 
in a normal direction from the surface. Beginning at the airfoil boundary, the 
derivatives Up: Vg are calculated by one-sided differences, and the metrics and Jacobian 
are again computed as described previouslv. These values are all local to subroutine 
EDDY. Ihe no-slip condition is now applied to determine the viscous stress, t, and 
skin friction. The friction value, CF, is one of the selectable normal output values 
from MAIN. The Van Driest damping factor, mr. is calculated and the 
subroutine begins scanning the grid points outward from the surface. The values ue, 
VE» Un Vu? the metrics and Jacobian; and Cae y, and F are calculated, and the values 
ol and v 


max "max 
the user-specified value ALFAI and 15 increasing or constant, the modified turbulence 


are found (Eqns. 2.27 and 2.28). If the angle of attack is less than 


model is used and F pass 1s the value of F where F* v is maximum. Above ALFAI on 
the upstroke and on the entire downstroke, the original Baldwin-Lomax model is used 
and D is the maximum value of F. Then the length scale, €, is found using the Van 
Driest function, and the inner-layer viscosity 1s computed (Eqn. 2.24). Next, the outer- 
laver formula previously described is used to find a value for viscosity at each normal 
grid location; and the point where the inner- and outer-layer values first cross 1s found. 
The final step performed by EDDY is to add the laminar viscosity to the appropriate 
inner- or outer-laver turbulent viscositv. 

Now that the viscosity has been calculated, subroutine STRESS can carry out 
evaluation of (Og R* + 0,5%) at the known time level, where R* and S* are, as given 
bv Equation 2.12, combinations of the vector R and S. First the derivatives in the 6- 
direction are computed by central differences. A more efficient discretization of the 
Spatial derivatives of viscous terms is realized by evaluation at half points, rather than 
at the nodal points. This reduces the number of nodal points involved from five to 
three. [Ref. 7: p. 26] Therefore, STRESS computes derivatives, metrics, Jacobian, and 


an average viscosity at the half points. These values are used to determine т.., US 
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yr Ry 


dissipation) and conduction (diffusion) terms, are set equal to zero. Finally the 


and S,. In the present implementation R, and S,, the viscous work (energy 


differenced value of де R* is added to the DQ1-DQ4 arrays. The same steps are then 
followed for the n-derivative terms. 

Subroutine SLPS has now assembled in the arrays DQ1-DQ4, the elements of the 
right-hand side of Equation 2.20 given by Equation 2.18, with explicit damping added 
but without the factor At. Next, the subroutine AMATI is called to compute the 
Jacobian matrix [A] = Golf" and the left-hand side of Equation 2.20 is then assembled. 
The second-order implicit damping terms are computed and added to the left side. 
Finally, the right-hand side (DQ arrays) is multiplied bv At and stored in matrix GG, 
giving the final form of the mght-hand side of Equation 2.20. A standard matrix 
inversion routine, MATRXI, obtains the solution to Equation 2.20, and Step 1 of the 
Beam-Warming algorithm, the $-sweep, is now complete. The solution, corresponding 
о Аф“ in Equation 2.20, is now stored in ће DQI-DQ4 arrays to become the right- 
hand side of Equation 2.21. The steps used in the -sweep are now repeated for the y- 
direction. AMAT2 computes the Jacobian matrix [B] = 0 G*, and the left-hand side 
of Equation 2.21 is assembled with implicit damping again added. MATRX2 1s called 
to obtain the solution (Aq*)n *! and complete Step 2 of the Beam-Warming algorithm. 
Step 3 is accomplished by updating the flow variables (Eqn. 2.22). Values at the 
boundaries are not updated. The outer boundaries remain at undisturbed freestream 
values making a specific step to apply boundary conditions unnecessary. This 
completes the solution procedure, but as a monitoring feature, the density elements of 
Aq* (now stored in DQ1) are scanned for the largest value. This value, along with its 
grid coordinates, 1s output at specified intervals, giving an indication of the most 
rapidly changing area of the flow field and the magnitude of the density change there. 

A solution has now been obtained by program MAIN. It remains only to apply 
the surface boundary conditions. This is accomplished by subroutine WALLBC. 
Boundary conditions were discussed in the last chapter, and the implementation is 
Straightforward. Conditions along the cut are averaged from above and below. The 
contravariant velocity component U at the two points nearest to the surface is used to 
determine U at the surface for Euler calculations; for viscous flows U = 0 at the 
surface. The appearance of the calculation steps is complicated slightly by use of a 
scheme that is applicable to both viscous and inviscid flows. Surface flow properties 


are finally obtained by extrapolating from interior points. 
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IV. RESULTS AND DISCUSSION 


Nine simulations are described in this chapter. Several variations were tried to 


investigate the effects of varying input parameters, some of which are mentioned 


without detailed presentation of results. The intention. was to obtain complete 


solutions under a range of conditions of dynamic effects, Reynolds number, Mach 


number, and turbulence model, whether or not such solutions are the most accurate 


obtainable under those conditions. [t was hoped thereby to gain some understanding 


of the most effective employment of the code, its limitations, and possibilities for 


improving its performance. Table 4 lists the input parameters for the simulations. The 


calculations were grouped into the following areas: 


1. 


2 
- 


(,2 


4. 


Steadv-state calculations at and bevond maximum lift. 


Comparison with deep stall experimental data, using the original Baldwin- 
Lomax turbulence model and modification. 


Comparison to experimental data for fully attached flow using both turbulence 
model versions. 


Low Reynolds number predictions at Mach numbers of .3 and .5. 


All of the runs were made with the implicit damping coefficient £j at a value of twenty 


times the explicit value £p (input as "WW"). The distance of the first point off the 


solid boundary was set to 0.00005 for all but the steady-state calculations beyond 


maximum lift. 


Uy 


Three basic types of output were obtained: 


The flow-field solution values p, ІШІ D and e and the corresponding grid 
coordinates at twelve equal time intervals for each cycle, for use in the plotüng 
program PLOT3D bv Pieter Buning of NASA Ames. This routine was used tó 


produce contour plots of pressure, densitv, Mach number, and stream function. 


Surface pressure coefficients at 96 equal time intervals for each cycle. This data 
was used in the plotting routine CARPET, written by Rosalie Lefkowitz of 
Sterling Software, to produce pressure distribution carpet plots utilizing the 
DISSPLA graphics library. 


The normal text output file containing, for 96 points each cycle, the integrated 
lift, drag, and moment coetficients: pressure distributions: the values an 
locations of the largest solution increments after every ten time Мерз and skin 
friction values. These coefficients were plotted using the DISSPLA plotting 


package QPLOT available on the NASA Ames VAX computers. 


The use of equal time intervals for output clusters the data near maximum and 


minimum angles of attack, where interest is often greatest. The solution was also 


saved on Cray tapes, at twelve points each cycle, for future program restarts, so a more 
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thorough investigation can be conducted of any interval with minimal use of computer 


time. 


TABLE 4 
INPUT CONDITIONS FOR COMPUTER SIMULATIONS 


1. a = 12°, Mo = .3, Re = 1,000,000 


Baldwin-Lomax Turbulence Model 


а = 14°, М. = .3, Ве = 1,000,000 


Baldwin-Lomax Turbulence Model 


. 4.= 135%, Mo = .301, Re = 3,910,000 
Modified Turbulence Model 


= О 


15°, а = 10°, Мо = .283, Re = 3,450,000 


Baldwin-Lomax Turbulence Model 


= 15°, @, = 10°, Meo = .283, Re = 3,450,000 
Modified Turbulence Model 


= 8, a, = 10", M = 184, Re = 2,450,000 
Baldwin-Lomax Turbulence Model 


5 @ = 8°, a, = 10°, Мо = .184, Ве = 2,450,000 
Modified Turbulence Model 


Мар = 15", d, = 10%, Mag = .284, Re = 345,000 
Baldwin-Lomax Turbulence Model 


15°, @, = 10°, Мо = .5, Ве = 345,000 


Baldwin-Lomax Turbulence Model 
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A. STEADY-STATE CALCULATIONS 

Before beginning dynamic simulations, the performance of the program in 
constant angle of attack calculations was explored. A verification series of runs at 
increasing angles of attack was made using the conditions reported by Sankar in 
Reference 6. The results agree with Sankar's and are not reproduced here. ее 
notable discrepancy from experimental results is the underprediction of maximum lift 
using the Baldwin-Lomax turbulence model, with separation occurring between one 
and two degrees below experimental values. The local time-step option, incorporated 
since Sankar's results were published, produced converged solutions in fewer than 1500 
steps at even the higher angles of attack. In addition to the integrated coefficients and 
surface pressure distribution reported bv Sankar, the flow field was graphically 
analvzed. Figure 4.1 contains densitv, pressure, Mach number, and stream function 
plots at 12 degrees angle of attack. 

Next the behavior of the solution above static stall angles was investigated. The 
original Baldwin-Lomax model was used with a Mach number of .3 and a Reynolds 
number of one million at 14 degrees angle of attack. The time-accurate mode was 
invoked by using a reduced frequency of .002 with zero amplitude of oscillation. The 
step size was .005 and the distance to the first point off the wall was .0001. Figure 4.2 
presents contour plots of density and pressure coefficient after 5100, 5610, and 6120 
steps. Based on the Mach number and time step, this should represent a total 
movement of three chord lengths of a particle in the freestream. The vortex, however, 
remains on the airfoil and gradually diminishes in intensitv. The integrated lift 
coefficient steadily decreased throughout the run to a value of less than 0.6. 

The modified turbulence model was applied to the conditions of frame 4019 of 
Reference S (volume 3, page 46) at the experimentally observed maximum lift angle of 
attack of 13.5 degrees. The modified turbulence model provides excellent results. For 
a Revnolds number of 3.91 million and Mach number of .301, the local time-step 
option was used with an explicit viscosity input of WW = 2. After 1500 time steps the 
experimental maximum lift coefficient of 1.4 was obtained, and the peak pressure 
coefficient of 8.1 appears to be very close to the experimental value. Moment and drag 
coefficient values also seem reasonable. The theoretical pressure distribution at 
maximum lift is shown in Figure 4.3, compared to the experimental distribution taken 


from Reference 8. 
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Figure 4.1 Steadv-State Attached Flow 
M go 7 A Ве= 1.00 x 109, a — 12? 
Top--Density (1), Pressure (r) Bottom--Mach No. (1), Stream Function (r). 
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Figure 4.2 Steady-State Separated Flow 


= 14° 


3, Re= 1.00 x 10°, a 
Density (1) and Pressure (r) at 5100, 5610, and 6120 time steps. 
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Figure 4.3 Steady-State Pressure Distribution at Maximum Lift 
М = .301, Ве= 3.91 х 106, а= 13.5° 
Inset--Experimental Data (Ref. 8). 
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B. COMPARISON WITH DEEP STALL EXPERIMENTAL DATA 

The first dynamic case studied was under the experimental conditions of frame 
9218 of Reference 8 (volume 3, page 146). These are the same conditions as Sankar 
used for comparison in Reference 6. Sankar reported only the integrated lift, noment, 
and drag coefficients, however, and here the surface pressure distributions and the flow 
field about the airfoil were investigated. Present calculations also used the exact Mach 
number and full Revnolds number, and runs were made with both the original Baldwin- 
Lomax model and the modification. This experimental case closely corresponds to the 
model dynamic stall (Figure 1.1). showing the formation of a strong vortex, an increase 
in maximum lift, rapid variations in pitching moment, and the typical lift versus angle 
of attack hvsteresis loop. The program was run under the following conditions: 

Re = 3.45 x 10° 


Мо = .283 
_ ye? 
й = 13 
эс ° 
и, = 10 
k = 151 


For the first run the original Baldwin-Lomax turbulence model was used with an 
explicit damping input of 5. The time step was a constant .005 throughout, and data 
was saved beginning at the mean angle on the second oscillation cycle. Figure 4.4 
shows the lift and pitching moment coefficients plotted against angle of attack for both 
theoretical and experimental data. Figure 4.5, a plot of the pressure drag coefficient, is 
included for completeness; drag plots are not presented for the remaining cases, since 
thev provided no additional information for this study. Figure 4.6 contains a carpet 
plot of the surface pressure coefficients for theoretical data and a plot of experimental 
data taken from Reference 8. Figures 4.7 to 4.14 are flow-field plots made at twelve 
equal time intervals during the cycle. The results show excellent qualitative agreement 
with experimental data throughout the cvcle. The moment coefficient is consistently 
low during the upstroke, but reasonable quantitative agreement of the lift coefficient 1s 
obtained until near 18 degrees. Just as in the steady-state case, the Baldwin-Lomax 
model gives an early prediction of flow separation. The maximum lift does not reach 
as high a value as experimental results, since separation occurs at a lower angle of 
attack in the theoretical results; as is the case for the experimental data, the lift 
continues to increase after the pressure distribution, indicated in the carpet plot, begins 


to break down. At the beginning of the downstroke, as the vortex is shed olf the 


4S 


THEORETICAL 
А 


= = ъ < -- - - ж = = єъ gw 


2.9 


LIFT COEFFICIENT 





5.0 10.0 15.0 20.0 25.0 


ANGLE OF ATTACK 


MOMENT COEFFICIENT 





9.0 10.0 15.0 20.0 25.0 


ANGLE OF ATTACK 


Figure 4.4. Lift and Pitching Moment Coefficients--Baldwin-Lomax Model 
M 3o 7.283, Re 3.45 x 109, a 9 15?—10? cos(.3t). 
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Figure 4.5 Pressure Drag Coefficient--Baldwin-Lomax Model 
М о = .283, Ве= 3.45 х 100, а= 15°—10° cos(.3t). 


50 


-10.0 





© 
N Х/С 





x/C 


Figure 4.6 Surface Pressure Coefficient Carpet Plots--Baldwin-Lomax Model 
Mag 7.283, Re 3.45 X 109, a — 15*—10 ^ cos(.30). 


Top--theoretical, Bottom--experimental (Ref. 8). 
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Figure 4.7 Density Contour Plots, Upstroke--Baldwin-Lomax Model 
Mg 7283, Rez 3.45 105799 15 92009 COSS 
6.31, 9.95, 14.93, 19.99, 23.65, 25 degrees. 
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Figure 4.8 Density Contour Plots, Downstroke--Baldwin-Lomax Model 
M so = -283, Re = 3.45 x 109, a — 15?—10*? cos(.3t) 
ADOS OS [003 6.36, 5 degrees. 
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Figure 4.9 Pressure Contour Plots, Upstroke--Baldwin-Lomax Моде! 
Мь=.233. Ве= 3.45 х 10%, a= 15%—10%cos(.3t) 
6.31, 9.95, 14.93, 19.99, 23.65, 25 degrees. 
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Figure 4.10 Pressure Contour Plots, Downstroke--Baldwin- Lomax Model 
M ag 7.283, Re- 3.45 x 106, a = 15? —10? cos(.3t) 
23.67, 20.02, 15.03, 10.03, 6.36, 5 degrees. 
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Figure 4.11 Mach Number Contour Plots, Upstroke--Baldwin-Lomax Model 
M „о = .283, Re = 3.45 Х 10°, a= 15°10 °cos(.3t) 
6.31, 9.95, 14.93. 19.99 2365 05 lo 
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Figure 4.12 Mach Number Contour Plots, Downstroke--Baldwin-Lomax Model 
M 59 = 283, Re = 3.45 x 105, a 2 15?—10? cos(.3t) 
6702002 1505003 6.56, 5 degrees. 
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Figure 4.13 Stream Function Plots, Upstroke--Baldwin-Lomax Model 
M 3g 7.283, Re= 3.45 x 10°, a 9 15? —10" cos(.3t) 
6.31, 9.95, 14,93, 19:99,:23:65 2510 G9 EEES. 
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Stream Function Plots, Downstroke--Baldwin-Lomax Model 
M gg, 7.283, Ве= 3.45 х 1006, а = 15?—10? cos(.3t) 
23.67, 20.02, 15.03, 10.03, 6.36, 5 degrees. 
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trailing edge, there is a jump in all the experimental coefficients. This is not duplicated 
in the theoretical data. The flow-field plots reveal the reason for this difference, as the 
vortex is not shed immediately but remains on the trailing edge and is dissipated during 
the downstroke. The absence of strong leading-edge peaks in the carpet plot, until the 
upstroke, shows the effect of this discrepancy on the pressure distribution. The 
recovery from separation occurs gradually as the vortex dissipates. Skin friction values 
are negative over the upper surface until the downstroke begins, and full recovery is 
not observed until completion of the downstroke. 

A run was made under these same experimental conditions but using the 
modified turbulence model on the upstroke. It was intended to turn off the modified 
model above 19 degrees, as recommended bv Sankar, but due to an error in the code, it 
was used throughout the upstroke. On the first attempted run, execution terminated as 
the angle approached 19 degrees. The values of Aq* had begun growing rapidly near 
the leading edge, and when a negatve value was computed for the density, a math 
error occurred causing the computer program to terminate abnormallv. The surface 
pressure distribution had been showing unusually high leading-edge peaks, but the 
integrated coefficient values were close to experimental values when termination 
occurred. Reducing the time step in half to 0.0025 at 15 degrees allowed the solution 
to proceed to above 23 degrees, but execution was then terminated in the same 
manner. The calculations were restarted at this time using a time step of .005 and 
distance of the first point off the wall of .00005, but this time with the viscosity input, 
WW = 10. Reduced, but very strong leading-edge suction peaks still appeared, the 
integrated coefficients were not as large, and a completed solution was obtained 
without reducing the time step. Data was again saved after one and one-fourth 
oscillation cvcles. 

The run using the modified turbulence model shows marked differences from the 
run with the Baldwin-Lomax model. The lift, drag, and moment coefficients are 
plotted along with experimental values in Figure 4.15 Now the rapid drop in lift is 
delayed too long--until the maximum angle is passed. Maximum lift is still 
underpredicted, but the higher viscosity value might be expected to affect accuracy. 
The moment coefficient is again lower than experimental during the upstroke until the 
vortex movement that is not predicted by the code. Note, however, the agreement with 
experiment once the "theoretical" vortex has been generated. The carpet plot, Figure 


4.16, reveals interesting differences from Figure 4.6. Until past the mean angle of 15 
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degrees, the pressure distributions are nearly identical. In fact, the Baldwin-Lomax 
model has a slightly stronger suction peak (—7.7 to —7.5 for the modified model). After 
this point, however, the Baldwin-Lomax model shows the effect of the premature flow 
separation while the suction peaks for the modified model grow abnormally high. Note 
the different scales used in the experimental and theoretical carpet plots. Although the 
integrated lift coefficient begins to drop slightly nearing the maximum angle, the 
pressure distribution is still apparently undisturbed at 25 degrees, with some weakening 
of the suction peak. The situation changes abruptly as the downstroke begins and the 
modified turbulence model is replaced by the original Baldwin-Lomax model. The 
suction peak finally breaks down, and by 17.6 degrees on the downstroke, the 
distribution is again very similar to that with the Baldwin-Lomax model. During the 
upstroke, the first negative surface friction values appear near the trailing edge at 21.6 
degrees, and bv 22.9 degrees, the negative values extend over the upper surface. This is 
the only clear indication of the residual stresses that are so abruptly relieved when the 
airfoil pitches down. 

The flow-field plots are similar in many respects to those for the Baldwin-Lomax 
model except for the timing of the vortex formation. Figures 4.17 and 4.18 are flow- 
field plots at the maximum angle and at 23.67 degrees on the downstroke. At the 
maximum angle there is still no strong vortex formation, although a narrow 
recirculating region is indicated near the trailing edge. The rapid development and 
shedding of the vortex occurs during the downstroke. By 23.62 degrees, the vortex has 
progressed well along the upper surface. The sequence of events is reasonable, 
although the timing is not. 

The program was corrected and run again under the same experimental 
conditions with two approaches to the problem of maintaining stability at the higher 
angles. First, the solution was marched from the mean angle (using the solution saved 
in the previous run at 15 degrees) to the maximum angle with the time step left at a 
constant .005 and the viscosity increased to 10. Then the run was repeated with the 
viscosity input left at 5 and the time step decreased to .0025 (as had been attempted 
earlier with the faulty turbulence model) This latter run was continued into the 
downstroke using a time step of .005, and the integrated lift and pitching moment 
coefficients are plotted in Figure 4.19. The results are improved over those with the 
Baldwin-Lomax model (Fig. 4.4). The maximum lift of 2.1 matches experiment and the 


lift curve slope is sustained longer into the upstroke. The negative moment peak is 
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Figure 4.15 Lift and Pitching Moment Coefficients--Modified Model until 258 
M 4g 2.283, Rez 3.45 X 109, a — 15?—10? cos(.3t). 


62 


-20.0 


-Ср 





2.0 


Х/С 





x/C 


с . - 0 
Figure 4.16 Surface Pressure Coefficient Carpet Plots-- Modified Model until 25 
Mo = .283, Re= 3.45 x 10°, a= | 8102-02-21) 


Top--theoretical, Bottom--experimental. 
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Figure 4.17 Density (top) and Pressure (bottom) Plots-- Modified Model until 258 
Мо=.283, Ве= 3.45 х 10%, а= 15?—10^ cos(.3t) 
25941) 25620850 
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Figure 4.18 Mach Number and Stream Function Plots--Modified Model until D 
M 4g = .283, Re = 3.15 x 109, a= 15*—10* cos(.3t) 
D T 
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Figure 4.19 Lift and Pitching Moment Coefficients-- Modified Model until 19° 
Mag 7.283, Re= 3.45 x 106, a = 15?—10* cos(.3t). 
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sharper. With the time step set at .005 and higher damping (WW = 10), results were 
not as good guantitatively. Peak values were reached at near the same points, but lift 
reached only 1.94, and peak moment was badlv underpredicted at —2.5. For a final 
comparison, the solution was restarted again for a short run from 15 degrees, but this 
time with the modified turbulence model turned off immediately. In all cases, a second 
leading-edge suction peak,indicative of vortex formation, was evident in the first output 
after the model was switched off. It seems from these results that the use of a smaller 
time step during the high angle portion of the upstroke is worth the added 


computational expense. 


C. COMPARISON WITH ATTACHED FLOW EXPERIMENTAL DATA 
The program was next run under the experimental conditions of Frame 9118 of 
Reference S. These conditions were chosen because, in contrast to the previous 
experimental results, flow remained attached throughout the oscillation cycle. 
Sufficiently lowering the reduced frequency (as in Frame 9110 of Reference 8) would 
result in separated flow. The program was run under the following conditions: 
Re = 2.45 x 10° 


E - га 
үзе о 

Uy = 8 

и, = 10° 

К = .199 


Again both turbulence models were tried, but for these conditions it was possible to 
complete both runs with a time step of .005 and explicit viscositv coefficient of 3. 
Output was again obtained after one and one-fourth cycles. 

As was expected, the Baldwin-Lomax model gave a separated flow prediction 
during the upstroke. The first negative surface friction values appear near the trailing 
edge at 16.3 degrees, and near 17 degrees a vortex forms at the leading edge. The lift 
and moment plots of Figure 4.20 show this discrepancy from experimental results. The 
lift falls below experimental values at high angles, but again, the action of the vortex 
sustains the lift coefficient and even causes the slope of the lift curve to increase 
slightly before the downstroke begins. The downstroke shows similar behavior to the 
dvnanmic stall case investigated previously, with a sharp drop in lift, which remains at a 
low value until the upstroke begins. The moment coefficient is again low during the 


upstroke and oscillates as the vortex propagates along the airfoil surface. The carpet 


67 


THEORETICAL 


LIFT COEFFICIENT 





-5.0 0.0 5.0 10.0 15.0 20.0 


ANGLE OF ATTACK 


0.05 


MOMENT COEFFICIENT 


20:15 





-0:20 
-5.0 0.0 5.0 10.0 15.0 20.0 


ANGLE OF ATTACK 
Figure 4.20 Lift and Pitching Moment Coefficients--Baldwin-Lomax Model 
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Mg 7.184, Re2 2.45 x 109, a= 8*—10*? cos(.4t). 


plot, Figure 4.21, shows the familiar ripple pattern of vortex movement. The flow field 
at the top of the cycle (18 degrees) is illustrated in Figure 4.22. 

Based on the previous results with the modified turbulence model, it was expected 
to perform better under this set of conditions than the Baldwin-Lomax model. The 
flow does in fact remain attached during the upstroke but it separates as soon as the 
downstroke begins, as it did in the dynamic stall case when the modified model was 
used throughout the upstroke. The lift and drag curves of Figure 4.23 are very similar 
in appearance to those of Figure 4.20 except for the loops formed at the maximum 
angle by the rapid vortex formation on the downstroke. The carpet plot, Figure 4.24, 
agrees quite well with experiment until the maximum angle. The flow field 1s shown by 
Figure 4.25 to be very similar at 16.64 degrees on the downstroke to that for the 


Baldwin-Lomax model at the maximum angle (Fig. 4.22). 


D. SIMULATIONS UNDER PROPOSED EXPERIMENTAL CONDITIONS 

A series of dynamic stall wind tunnel experiments 1s presently being planned at 
the Ames Research Center Fluid Mechanics Laboratory. These experiments, which 
will include investigation of compressibility effects, are to be conducted in a Revnolds 
number range lower by a factor of ten than those reported in Reference 8. The last 
two simulations attempted, therefore, were made at a Reynolds number of 345,000. In 
order to have some basis for comparison, the mean angle, amplitude of oscillation. and 
reduced frequency chosen were the same as for the deep dynamic stall, high Reynolds 


number runs. The Mach numbers used were .284 and .5. The conditions then were the 


following: 
Re = 345,000 
М. = .284 (Run 1), .5 (Run 2) 
ар = 15° 
а, = 10° 
К = .151 


The conditions are within the range of the planned experiments. The Baldwin-Lomax 
model was used for both runs. The first run was conducted with a constant time step 
of .005 and viscosity coefficient input of 5. At the higher Mach number of .5, the time 
step was set to .0025. An attempt to start the run with a .005 time step resulted in 


abnormal termination after fewer than 30 steps due to the calculation of negative 


69 


% deg. 
“Л ‚68 
-0.69 


“20.0 





220 


ah” 
O 
O 


Х/С 





x / C 


Figure 4.21 Surface Pressure Coefficient Carpet Plot--Baldwin-Lomax Model 
М о = .184, Ке= 2.45 х 10, а = 8° 10° соѕ(.40) 


Top--theoretical, Bottom--experimental (Ref. 8). 
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Figure 4.22 Flow-field Plots-- 
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Figure 4.23 
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Lift and Pitching Moment Coefficients-- Modified Turbulence Model 
Му = .184, Ве= 2.45 х 100, а= 8°—10°соѕ(.40). 
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Figure 4.24 Surface Pressure Coefficient Carpet Plots--Modified Turbulence Model 
Mo = .184, Re= 2.45 x 105, a= 8°—10°cos(.dt) 


Top--theoretical,Bottom--experimental (Ref. 8). 
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Figure 4.25 Flow-field Plots-- Modified Turbulence Model 
М. 7.184, Rez 2.45 x 109, a — 8*—10? cos(.-It) 
lop--Density (1), Pressure (r) Bottom-- Mach No. (lI), Stream Function (r). 
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density near the leading edge, just as had happened during the earlier high Reynolds 
number run with the modified turbulence model. With the smaller time step, little 
improvement was realized until the artificial viscosity input was increased. Fewer than 
fifty steps were completed with the explicit viscosity at 5. When this was doubled the 
solution reached nearlv 8 degrees angle of attack and completed 2170 steps before it 
again “blew up”. In each instance in which this abnormal termination took place, the 
output values all appeared normal until the values of the residuals began growing 
shortly before ternunation. Increasing the explicit viscosity coefficient to 20 permitted 
a complete solution to be obtained. 

The results for these two runs were qualitatively similar to those for the higher 
Reynolds number. The lift and moment coefficients plotted in Figure 4.26 show values 
close to those of the earlier results. The slopes of the lift curves are somewhat steeper 
than the previous results, and the lift continues to drop (and the moment coefficient to 
rise) farther into the downstroke. Figure 4.27 compares the pressure distribution at the 
lower Mach number with earlier results for the high Mach number. The suction peaks 
at the low Reynolds number are not as strong and begin breaking down slightly earlier. 
The flow-field plots showed a sinular series of events to those seen previously and are 
motpresented here. 

The results obtained at the higher Mach number have two notable features: 
smoother contours and lower peak values. At least some of this smoothing must be 
attributed to the much higher viscosity input used for this case (20 versus 5). The lift 
and moment plots are given in Figure 4.28 and the carpet plot for this case along with 
that at the lower Mach number are shown for comparison in Figure 4.29. The shape 
of Figure 4.28 curves is more rounded than was Figure 4.26, and maximum lift is 
reduced. The pressure coefficient suction peaks (Fig. 4.29) are also smoother, lower, 
and begin breaking down sooner at the higher Mach number. To further investigate 
the effect of artificial viscosity, a partial run was made with the explicit damping 
coeflicient input as 55. The maximum lift coefficient obtained was 1.3 at 17.6 degrees 
angle of attack compared to 1.6 at 18.2 degrees with the viscosity at 20. The maximum 
value was attained at nearly the same angle (output was not taken at exactly the same 
points for this second run). However, the first irregularities in the leading-edge 
pressure coefficient were delaved by this increase in viscosity from about 12.3 degrees 
Eo attack to near 15 degrees. Thus, the computer results do seem to be 
predicting an earlier vortex formation (and breakdown in the surface pressure 


coefficient) at the higher Mach number. 
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Figure 4.26 Lift and Pitching Moment Coefficients--Baldwin-Lomax Model 
M9 = .284, Re=.345 x 10°, a= 15°-10° cos(.3t). 
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Figure 4.27 Surface Pressure Coefficient Carpet Plots--Baldwin-Lomax Model 
Top: May =.284, ЕВе-.345 х 106, а-157-107со5(.30 
Bottom: Re= 3.45 X 10°(Previous Results). 
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Figure 4.28 Lift and Pitching Moment Coefficients--Baldwin-Lomax Model 
M œ =.5, Re=.345 x 106 a= D W 8! 
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Figure 4.29 Surface Pressure Coefficient Carpet Plots--Baldwin-Lomax Model 
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V. CONCLUDING REMARKS 


The computer program emploved for this study has shown the capability to 
model the basic events of the dynamic stall process. The quality of the results is, not 
surprisingly, highest at moderate angles--two to three degrees above static stall angle in 
examples studied. Results at higher angles are strongly dependent on turbulent 
viscosity model. The downstroke is not as well approximated, although there is 
qualitative agreement of the integrated coefficients. The program appears quite robust 
at Mach numbers of .3 and below using the Baldwin-Lomax turbulence model. With 
the modification to the turbulence model introduced in this implementation, better 
results are obtainable at the higher angles, but the program is less stable. Adequate 
convergence was generally achieved quickly enough that data could be saved beginning 
at the mean angle on the first oscillation cycle. 

The examples considered show some of the limitations of current Navier-Stokes 
solvers. By manipulating artificial viscosity and turbulence modelling. it is possible to 
manage calculations to match, in some degree, desired results. For example, the 
attached flow case considered might be brought into closer agreement with experiment 
by using the modified turbulence model after the beginning of the downstroke. For the 
dynamic stall case, however, this might lead to premature flow reattaclunent and less 
accurate results. The simple modification introduced, used as recommended during the 
upstroke, shows encouraging results. At no significant expense, modelling near and a 
few degrees above static stall angles appears improved without affecting results at 
lower angles. Further testing needs to be done to determine the effectiveness of the 
change under different dynamic conditions, and the use of other models is certainly 
worth consideration. King and Johnson have reported favorable results with a 
nonequilibrium turbulence model using an ordinary differential equation to model 
streamwise shear stress development, and an eddy viscosity distribution across the 
boundary laver based on the maximum value from this equation. Compared to the 
Baldwin-Lomax and Cebeci-Smith models, it showed little difference in attached flow 
calculations with a Navier-Stokes solver using a Beam-Warming scheme. When 
applied to separated flows and flows with strong viscous effects, it gave superior 


results, though at some computational expense in its present form. [Refs. 18,20] 
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The use of other grids or variations on the present one would be worth 
considering for the purpose of developing comparisons for the proposed experiments at 
the Ames Research Center Fluid Mechanics Laboratory. Obtaining a solution at the .5 
Mach number was difficult and required a large increase in the artificial viscosity. 
Some of the problems might be alleviated by closer spacing in the leading-edge area 
where the instability arose. The present grid generation scheme might be tried with 
more points selected in the normal direction. Events at the trailing edge are not now 
modelled with complete accuracy. The vortex moves to the trailing edge and diffuses 
there. This behavior may be influenced by the level of artificial viscosity (Which must 
be increased to compensate for grid coarseness), the downstream boundary, location of 
the grid cut, and relative coarseness of the grid at the trailing edge. The present grid 
provides adequate results for a large range of applications with a reasonable number of 
grid points, however. 

The solver now assumes a fully turbulent boundary laver (or fullv laminar if the 
Reynolds number is set to zero). The proposed experiments may include significant 
transition point effects. The transition could be simulated by retaining the laminar 
viscosity. coefficient forward of a specified chord location. This might be especially 
useful for comparison with tripped boundary layer data. Some runs were made with 
fully laminar viscositv. The pressure gradients at the leading edge could not be 
sustained, and multiple vortices were formed but the solution remained stable. 

The cases considered in this study represent useful test conditions for studying 
the effect of changing various program segments. The deep stall case has the model 
dynamic stall features and has been used by Sankar. The attached flow case is an 
interesting and a difficult contrasting case. The forthcoming experiments will provide 
the opportunity to investigate details of the dynamic stall process, including 
compressibility effects, in much greater depth. If the Navier-Stokes solver is not yet a 
completely adequate predictive tool, the effort to make it so can, of itself, provide 


insight into the nature of unsteady flows. 
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APPENDIX A 
SANKAR NAVIER-STOKES SOLVER 


JOB, JN=, T=. 
ACCOUNT , AC=, US=, UPW=. 


.. 
СҒТ,ОМ-А,ОҒҒ-5. 


°. 

. ASSIGN , DN=NEWS LN , A=F T08. 

. ASS IGN , DN=NEWCP , A=F T18. 

. ASSIGN , DN=XYZ , A=FT11. 

. ASSIGN, DN=Q, A=FT12. 

.IS THIS RUN TO START FROM A STORED SOLUTION? 

. ACCESS , DN=OLDSLN , PONm , I D= , 

.DO YOU WANT TO ACCESS EXISTING PRESSURE COEFFICIENT DATA? 
.ACCESS , DN=OLDCP , PDN=, ID=. 


"eo 9 0 909 @ @ #é еее €. 9 9 9 9 € e € e e 9 ег ет 9? 9 е е е еее e e 


ASSIGN, DN=OLDSLN, A=FTQ7. 
ASSIGN, DN=OLDCP, A=FT17. 


DR ,MAP=PART ,SET#=ZERO. 

.DO YOU WANT TO SAVE THE CURRENT SOLUTION? 

. SAVE , DN=NEWSLN, PON=, ID=. 

.DO YOU WANT TO SAVE PRESSURE COEFFICIENT DATA? 
. (PREVIOUS FILES SHOULD BE PURGED FREQUENTLY). 
. SAVE , DN=NEWCP , PDN=, [D=. 


DO YOU WANT TO CREATE PLOTSD FILES? 

. ACCESS , DN=SENDVAX , PON=SENDVAX , IDzSTTRDM. 
. SENDVAX , DN=XYZ, VDN=. 

. SENDVAX , DN=Q, VDN=. 


.DO YOU WANT PLOT3D FILES FOR TROUBLESHOOTING IF PROGRAM CRASHES? 
EXIT. 

. ACCESS , DN=SENDVAX , PON=SENDVAX , ID=STTRDM. 

. SENDVAX , ON=XYZ,VDN=. 

. SENDVAX , DN=0Q , VDN=. 


.DO YOU WANT TO USE JOB CHAINING? 
. РЕТСН, Омө, ТЕХТэ. 

. REWIND, DN». 

. SUBMIT, ON=. 

EOF 


PROGRAM MAIN 

PROGRAM MAIN( INPUT, OUTPUT, TAPES=INPUT , TAPE6=OUTPUT ) 
COMMON /SURF/PSUR( 161) 

COMMON / F IX /OMEGA 

COMMON /MUTUR/CMU( 161,41) 

COMMON/SKINCF/CF(161) 
COMMON/GRID1/X(161,41),2(161,41) 
COMMON/PAR/GAMMA , REYREF , ALFA, ALFA1,REDFRE, AMINF , ALFAIT 
COMMON/DGRID/DT , IMAX,KMAX, ITEL, ITEU 


COMMON/GR ID/YACOB( 161,41) 
COMMON/DAMP /WW , WW T 
COMMON/FLOW/Q1(161,41),Q2(161,41),Q3(161,41),04(161,41) 
DIMENSION TITLE(15), TYTLE(15) , ALPHA(96) ,CPTH(97,96) , XTH(97) 
COMMON/LOGIC/RSTRT,PITCH,PLUNGE 
LOGICAL RSTRT.PITCH,PLUNGE 
С 
C!!!!!INDICATES COMMENTS OR CHANGES ADDED BY J.F. VALDES, DEC 1986 
Ceee PROGRAM SOLVES TWO-DIMENSIONAL FLOW PAST ARBITRARY 
Cees GEOMETRIES USING ADI PROCEDURE. 


TAPES = FILE CONTAINING INPUT DATA 

TAPE6 = OUTPUT 

TAPES = FILE THAT SAVES THE FLOW FIELD AT THE END OF A RUN 

IF THE CURRENT RUN IS A RESTART OF A PREVIOUS RUN, THEN 

TAPE7 IS USED TO READ THE FLOW FIELD INTO MEMORY 

ITAPE18 IS USED TO ACCUMULATE PRESSURE COEFFICIENT DATA FOR A CYCLE 
! IT MUST BE ACCESSED WHEN COMPLETE BY PROGRAM PLOTNSE.SRC 

I TAPE17 IS USED TO READ EXISTING PRESSURE COEFFICIENT DATA FOR 
! UPDATE DURING THE CURRENT PROGRAM RUN 
I TAPE11 [S USED FOR PLOT3D XYZ-FILE OUTPUT 
! TAPE12 IS USED FOR PLOT3D Q-FILE OUTPUT 


geff w geg += gem = 
cw w gem cum ge фын» 


0000000000000 


Сә» READ INPUT DATA 
С 

PI = ATAN(1.)*4 

READ (5,1) TITLE 

READ (5,2) IMAX,KMAX,DT,WW,ALFA,ALFA1,ALFAI ,REDFRE,AMINF 


ALFAD = ALFA 


C 
C FORCE DT TO BE EQUAL TO UNITY FOR STEADY FLOW PROBLEMS 
C THIS INVOKES THE RELAXATION MODE 

C 


IF(REDFRE.LE.0.001) DT = 1.0 
READ (5,2001) 


C NSTP = NO. OF TIME STEPS TO BE DONE ON THIS RUN 
READ (5,2221) FNSTP 
NSTP = FNSTP 
READ (5,2001) 
С REYREF= REYNOLDS NUMBER IN MILLIONS 
С DNMIN = DISTANCE OF FIRST POINT OFF THE WALL 
С FOR REYNOLDS NUMBERS UPTO 3 MILLION USE 0.00005 
READ (5,2221) REYREF,DNMIN 
REYREF = REYREF + 1.E+06 
C 
C  TSTART = TIME THAT THE CALCULATIONS HAVE BEEN ADVANCED 
C UPTO THE PREVIOUS RUN. IF TSTART IS NEGATIVE THIS VALUE IS 
C X OBTAINED FROM TAPE 8. 
READ 22220 
READ (5,2221) TSTART 
CIIHFFIFULOUT IS A FLAG FOR GENERATING PLOTTING FILES. WHEN NEGATIVE 
C!!!!!NO DATA IS GENERATED. ZERO IS USED TO START AND A POSITIVE 
CI!!! INUMBER TO CONTINUE GENERATING FULL OUTPUT. ADDED FEATURE(JFV). 


READ 5. 2951) 
READ (5,2221) FULOUT 
2221 FORMAT(2F10.4) 
READ (5,2001) 
2001 FORMAT(1X) 
READ (5,2000) RSTRT,PITCH,PLUNGE 
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2000 FORMAT(4L5) 


C NEGATIVE REYREF MEANS INVISCID FLOW 
C 

Cese PRINT OUT THE INPUT DATA 

C 


WRITE (6,4) TITLE 
WRITE (6,3) IMAX,KMAX,DT,WW,ALFA, ALFAT,ALFAI,REDFRE, AMINF 
IF(REYREF.GT.0.) WRITE (6,3700) REYREF 


САММА=1 . 4 
C 
C!!!!!GENERATE COMPUTATIONAL GRID 
С 
CALL AIRFOL( IMAX,KMAX, ITEL, ITEU) 
IF (REYREF.GT.@. )CALL CLUSTR(DNMIN) 
WRITE (6,3002) 
C 


Cees STARTING CONDITIONS. 
Cees DENSITY NORMALISED WITH RESPECT TO ROINF 
Cees VELOCITIES NORMALISED WITH RESPECT TO AINF 
Cees TOTAL ENERGY NORMALISED WITH RESPECT TO (ROINFeAINFeAINF) 
C 
TOTENSAMINF *AMINF 0. 54-1. / (GAMMAs (GAMMA-1 . )) 
ALFA = АСРА « АТАМ(1.) / 45. 
ALMEAN = ALFA 
ALFAI = ALFAI * ATAN(1.) / 45. 
ALFA1 = ALFA1 в АТАМ(1.) / 45. 
ALFAS » ALMEAN - ALFA1 
UINF » AMINF e SN 
VINF = AMINF + SIN(ALFA) 
00 7 1«1,1МАХ 
DO 7 K=1,KMAX 
Q1(1,K)=1. 
Q2(1,K)=UINF 
Q3(1,K)=V INF 
Q4(1,K)=TOTEN 
7 CONTINUE 
IF(RSTRT) READ (7) TIME,Q1,Q2,03,Q4 
IF(TSTART.GE.@.) TIME = TSTART 
IF(.NOT.(RSTRT)) TIME = 0. 


C 
CI!1! H!READ STORED PRESSURE COEFFICIENTS. STATEMENTS ADDED BY JFV. 
C 
IF (FULOUT .GT. 0.) THEN 
READ (17,1) TYTLE 
READ (17) AMPLTD,STMEAN,OMS, XTH, ALPHA, CPTH 
END IF 
CI11! IBEGIN PRESSURE COEFFICIENT DATA FILE. STATEMENTS ADDED BY JFV. 
IF (FULOUT .EQ. 0.) THEN 
АМРЕТО » ALFA1 * 45. / АТАМ(1.) 
STMEAN = ALFA • 45. / ATAN(1.) 
OMS = REDFRE 
TEDGE = X(ITEL+48,1) 
DO 15 I=1,97 
XTH(I) = X(I+ITEL-1,1) - TEDGE 
15 CONTINUE 
END IF 
Cl!!! 1DETERMINE NUMBER OF TIME STEPS OF CURRENT SIZE IN ONE CYCLE 
Citii! SCHEME FOR PLOTTING DATA AND PROGRAM EXIT ADDED BY ЈЕУ 
IF (PITCH) CYCLE = PI / (REDFRE « AMINF » DT) 
C 224 
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о 


CALL METRIC 
DO 1000 ITN=1,NSTP 


C 
CIIIIICREATE PLOT3D FILES FOR INVESTIGATING INSTABILITY ADDED BY JFV 
СІРІ! THIS 8LOCK UPDATES PLOT3D OUTPUT FILE AFTER 
СІРІ EACH STEP. USE IN CONJUNCTION WITH EXIT JCL. 
С11!1! sss» NOTE — USE ONLY FOR TROUBLESHOOTING eee 
C 
C REWIND(11) 
C WRITE(11) IMAX, KMAX 
C WRITE(11) ((X(I,J), I21,IMAX), J=1,KMAX), 
C 1 ((Z(I,J), I=1,IMAX), J=1,KMAX) 
C REWIND( 12) 
С WriTE(12) ЇМАХ, КМАХ 
С WRITE 12) AMINF, ALFAD, REYREF, TIME 
С WRITE(12 E I,J), 1=1,1МАХ), /=1,КМАХ), 
С 1 шиг) 141,1МАХ), /=1,КМАХ), 
С 2 ((Q3(1,J), I=1,IMAX), J=1,KMAX), 
C 3 ((Q4(I,J), I21,IMAX), J=1,KMAX) 
C 
TIME = TIME + DT 
CIIITIROTATE GRID TO NEW ANGLE OR SET TO CORRECT ANGLE FOR RESTARTS 


IF (PITCH) THEN 
OMEGA = 2. » REDFRE *AMINFeSIN(REDFRE » 2. » TIME * AMINF) 


1 *ALFA1 
ALOLD = ALMEAN — ALFA1 » С05(2. + REDFRE + AMINF e 
1 (TIME — DT)) 


ALFA = ALMEAN — ALFA1 + COS(REDFRE + 2. + TIME * AMINF) 
ALFAD = ALFA e 45. / ATAN(1.) 
DALFA = ALFA — ALOLD 
IF(RSTRT.AND.ITN.EQ.1) DALFA = ALFA — ( ALMEAN - ALFA1 ) 
CALL ROTGRID(X,Z,IMAX,KMAX,DALFA) 
CALL METRIC 
END IF 
IF (PLUNGE) THEN 
OMEGA = 2. » REDFRE * AMINF 
ALFA - А(МЕАМ 
END IF 


C 
C!!!!!COMPUTE THE SOLUTION BY ADI TECHNIQUE 
CALL SLPS(ITN,OMEGA,DALFA) 


APPLICATION OF BOUNDARY CONDITIONS 
CALL WALLBC 


C 
C 
C 
C 
C PRINT OUT PRESSURE AT THE SURFACE 
С 
С!!!!!ЕОВ $ТЕАОҮ $ТАТЕ ООТРОТ USE THE FOLLOWING 
IF((1TN/10@*100.EQ.1TN) .AND. (.NOT.PITCH)) THEN 
WRITE юэ 
WRITE (6,33) ITN 
CALL CPPLOT(ITEL, ITEU, AMINF,X(1, 1), Z( 1, 1) , PSUR) 
WRITE (6,12) (I,CF(I), I», IMAX) 
CALL LOAD(PSUR,CL CD, CM, ALFAS) 
WRITE (6,3000) CL, CD, CM 
WRITE (6,3002) 
END IF 
C!!11!0R, FOR DYNAMIC SOLUTION OUTPUT AT EQUAL PHASE INTERVALS 
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C!!!!! MODIFIED OUTPUT SCHEME AND PROGRAM EXIT SCHEME BY JFV. 
IF (PITCH) THEN 
TINCR = TIME/OT - CYCLE/4. 
IF (TINCR .LT. 0.) TINCR » TINCR + CYCLE 
CI!!! DETERMINE NUMBER OF STEPS BETWEEN OUTPUT 
NCPOUT = NINT(CYCLE/96. ) 
CI!!! !MULTIPLY NCPOUT BY THE NUMBER OF OUTPUT CYCLES DESIRED BETWEEN 
C!!!!!PROGRAM EXITS. THIS DETERMINES THE INTERVAL FOR PLOT3D FILES. 
NEXIT = 24 * NCPOUT 
ACYC = AMOO(TINCR,CYCLE) 


DO 19 J=1,4 

IF (ABS(TINCR - CYCLE+J) .LT. 0.5 .AND. (ITN .GT. 

1 NCPOUT)) ACYC = 96. / (NEXIT / NCPOUT) + NEXIT 
19 CONTINUE 


NBCYC = NINT(ACYC) 

IF (NBCYC/NCPOUT + NCPOUT .EQ. NBCYC) THEN 
INDEX = NBCYC / NCPOUT 
ALPHA( INDEX) = ALFAD 


C .... 
C 
WRITE (ШІ 
WRITE (6,33) ITN 
WRITE (6,3500) ALFAD 
CALL CPPLOT(ITEL, ITEU, AMINF,X(1,1),2(1, 1) , PSUR) 
C 
CIII!I!ISTORE CURRENT PRESSURE COEFFICIENTS. ADDED STATEMENTS BY JFV. 


DO 20 J=1,97 
CPTH(J, INDEX) = PSUR(J+ITEL-1) 
20 CONTINUE 


C 
СЕМ 
IF (FULOUT .GE. 0.) WRITE(6,12) (I,CF(I),Is1, IMAX) 
CALL LOAD(PSUR,CL, CD, CM, ALFAS) 
WRITE (6,3000) CL , CD , CM 
C 
C!!!!!FOR AUTOMATIC PROGRAM EXIT DURING FINAL OUTPUT (JFV) 
C 
IF ((NBCYC/NEXIT + NEXIT .EQ. NBCYC) .AND. (ITN .GT. 
1 2 , МСРООТ)) СО ТО 1901 
WRITE (6,3002) 
C 994. 
С 
END IF 
END IF 
C 
C 
1900 CONTINUE 
СРЇН 


1991 CONTINUE 
WRITE (8) Т1МЕ,01,02,03,04 
CIII!!SAVE PRESSURE COEFFICIENT DATA (JFV). 
IF (FULOUT .GE. 0.) THEN 
WRITE A TITLE 
WRITE (18) AMPLTD,STMEAN, OMS, XTH, ALPHA, CPTH 


I! 1T ICREATE PLOTSO FILES. FEATURE ADDED BY УРУ. 
11 1 11DO NOT USE WHEN TROUBLESHOOTING (ABOVE) 


ОООО 


WRITE(11) IMAX, KMAX 
WRITE(11) ((X(I,J), I21, IMAX), J=1,KMAX), 
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1 ((2(1,4), 141,1МАХ), 4а!,КМАХ) 
WRITE(12) IMAX, KMAX 
WRITE(12) AMINF, ALFAD, REYREF, TIME 
WRITE(12) 1118! 141, МАХ), NAT. 
1 ((02(1,4), 1=1,1МАХ), 4=1,КМАХ), 
2 ((03(1,4), 1=1,1МАХ), in E 
((04(1,4). 1=1,1МАХ), J=1,KMAX 


ЕМО ЇЕ 
С 
С PRINT OUT VELOCITY PROFILE 
C 
DO 4000 I = 1 , IMAX , 2 
5 20. 
DO 4000 K 22 , 10 
S = S + SORT((X(I,K)-X(I,K-1)) **2*(Z(I,K)-2(I,K-1))*92) 
ED = CMU(I,K) / DT 
UTOT » SQRT(Q2(I,K)*»2 «4 03(1,К)9д2)/01(1,К) 
C!!!!!IF PRINTED OUTPUT IS DESIRED 
C WRITE (6,2002) I , K,Q2(I,K) , Q3(I,K) , ED , UTOT,S 
4000 CONTINUE 
STOP 
1 mA 315 
2 FORMAT(215,7F10.0) 
3 FORMAT(/,5X, SHIMAX=, 110, //,5X, SHKMAX=,110,/,5X,5H OT=,F20.8, 


1/,5X,5H W=,F20.8,/,5X,SHALFA=,F20.8,/,5X,6HALFA1=,F20.8, 
2/,5X, 6HALFAI=,F20.8,/,5X, 7HREDFRE=,F20.8,/,5X, SHAMINF=,F20.8) 
4 RCRMAT (171 8X: 154) 
12 FORMAT(8(14,F10.4)) 
19 FORMAT(/) 
22 FORMAT(2F10.6,15) 
33 FORMAT(5X, 5HISTP=, I5) 
2002 FORMAT(215,5E14.6) 
3000 FORMAT(5X,3HCLs,F10.4,5X, 3HCDs,F 10.4, 5X, 3HCMs,F 10.4) 
3002 FORMAT(//,4X, 'DRMAX' , 11X, 'DUMAX' , 11X, 'DVMAX' , 11X, 'DEMAX' , 1OX, 
1 IR^ ,3X, 'KR') 
3500 A оа) 
3700 FORMAT(5X,7HREYREF=,F13.4) 
END 
SUBROUTINE AMAT1(K,IMX1,XIX,XIZ,XIT) 
COMMON/F LOW/Q1(161,41),02(161,41),03(161,41),04(161, 41) 
COMMON/PERTR/DQ1 (161,41),002(161,41),003(161,41),004(161,41) 
COMMON/AM/A(4, 4,161) 
COMMON/PAR/GAMMA , REYREF , ALFA, ALFAT, REDFRE, AMINF , ALFAI 
DIMENSION XIT(161,41),XIX(161,41),XIZ(161,41) 
REAL K@,K1,K2 
С 
Севе AMAT! COMPUTES THE COEFFICIENT MATRIX DE/DQ DURING XI SWEEP 
C 


GM1 = GAMMA — 1 

DO 1000 I = 2 , IMXI 

ке = E 

K1 = XIX(I,K) 

K2 = XIZ(I,K) 

U = 02(1,К) / Q1(I,K) 
W = ООШ, / Q1(1,K) 
ЕВУК = Q4(1,K) / ыс 
РНІ2 = 0.5 CM1 $» (U » U+ We W) 
THETA = К1 * U + K2 * W 
Nor] - KO 

А(1,2,1 = K1 
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1000 


Сеге 


1999 


КІ « PHI2 — U « THETA 

КӨ + ТНЕТА - К! * (61 -1.) , 0 

K2 * U - ОМ1 s KI e W 

GM1 

РН12 – W « THETA 

W- K2 •« СМІ e. U 

KO + THETA - K2 * (GMI - 1.) * W 

K2 в СМ1 

ТНЕТА в (2. * PHI2 - GAMMA * EBYR) 

Ki « (GAMMA • ЕВҮВ - PHI2) - GM1 * U * THETA 
K2 « (GAMMA « EBYR - PHI2) - GM! * W * THETA 
КӨ % САММА * THETA 


A 
N 
Ф о е 


H H H H H ll H H HMH M M M NM 
Ж 


ЕМО 

SUBROUTINE AMAT2(1,KMX1,ZETAX, ZETAZ, ZETAT ) 
COMMON/FLOW/Q1(161,41),Q02(161,41) ,Q3(161,41) ,04(161, 41) 
COMMON/PERTR/DQ1(161,41),DQ2(161,41),DQ3(161,41),DQ04(161,41) 
COMMON/AM/A(4,4,161) 

COMMON/PAR/GAMMA , REYREF , ALFA, ALFA1, REDFRE, AMINF , ALFAI 
DIMENSION ZETAX(161,41),ZETAZ(161,41),ZETAT(161,41) 

REAL K@,K1,K2 


AMAT2 COMPUTES THE COEFFICIENT MATRIX DF/DQ DURING ETA SWEEP 


GM1 z GAMMA - 1. 

DO 1000 K 2 2 , KMX1 

ке = ZETAT(I,K) 

K1 = ZETAX(I,K) 

K2 = ZETAZ(I,K) 

U = ЕШ! 1 10) 

W = Q3(I,K) / Q1(1,K) 

EBYR = Q4(1,K) / 01(1,K) 

PHI2 = 0.5 * GMI * (U * U + W * W) 

THETA = К! * U + K2 ° W 

A(1,1,K) = K@ 

A(1,2,K) = K1 

NAE = K2 

A 35 ж 0 

мас = Ki « PHI2 -U « THETA 

A(2,2,K) = K@ + THETA — K1 « (GM1-1.) • 0 
A(2,3,K) = К2 * U — GM1 e K1 «e W 

A(2,4,K) = K1 e GM1 

A(3,1,K) - K2 * PHI2 - W èe THETA 

A(3,2,K) = K1 s W- K2 e GM1 e U 

A(3,3,K) = КӨ + THETA -K2 * (GM1-1.) * W 
A(3,4,K) =з К2 • СМІ 

A(4,1,K) = THETA * (2. * PHI2 - GAMMA • ЕВҮВ) 
А(4,2,К) «КІ e (GAMMA « EBYR - PHI2) – СМ1 • U è THETA 
A(4,3,K) = К2 e (GAMMA « EBYR - PHI2) - GM1 * W e THETA 
A(4,4,K) = K@ + GAMMA « THETA 

CONTINUE 

RETURN 

END 


SUBROUTINE SLPS(ITN,OMEGA,DALFA) 

COMMON/F LOW/Q1(161,41),Q2(161,41),Q03(161,41),04(161,41) 
COMMON/PERTR/0Q1(161,41),002(161,41),003(161,41) ,004(161,41) 
СОММОМ /АМ/А( 4, 4, 161) 
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COMMON/TRID/DD(4,4,161,41),MM(4,4,161,41), EE(4,4,161,41) 
1,66(4,161,41) 
COMMON/PAR /GAMMA , REYREF , ALFA, ALFAT , REDFRE, AMINF , ALFAT 
COMMON/DGRID/DT , IMAX,KMAX, ITEL, ITEU 
COMMON/GRID/YACOB(161,41) 
COMMON/D AMP /WW , WWI 
COMMON/MTRIX/ ХІХ(161, 41), ХІ2(161, 41), 2ЕТАХ(161, 41), 2ЕТА2(161,41) 
1 ,XIT(161,41),ZETAT(161,41) 
REAL MM 
DIMENSION DELTAT(161,41) 
C 
Ceee SUBROUTINE SLPS DOES THE BULK OF THE WORK FOR THE ADI ALGORITHM. 
Ceee IT CALLS FLUX AND COMPUTES RIGHT HAND SIDE DURING THE TWO SWEEPS, 
Coro ASSEMBLES THE COEFFICIENT MARICES, ADDS IMPLICIT AND EXPLICIT 
Ceee DISSIPATION AND CALLS THE TRIDIAGONAL SOLVER TO OBTAIN THE FINAL 
Cesa SOLUTION. 
С111115ЕТ VALUE OF IMPLICIT DAMPING COEFFICIENT 
WWI = 20.0 o Wé 
IM1 = [МАХ - 1 
[М2 = [МАХ - 2 
КМ1 = КМАХ - 1 
КМ2 = КМАХ - 2 
IF(ITN.EQ.1) THEN 
CI!!!!LOCAL TIME STEP OPTION FOR STEADY STATE CALCULATIONS 
IF(REDFRE.LT.0.001) THEN 
DO 777 K = 2 , KMAX - 1 
ОО 777 1 2 , IMAX - 1 
DELTAT(I,K) 2 0.5 / (1. * SORT(ABS(YACOB(I,K)))) 
777 CONTINUE 
ELSE 
DO 778 K = 2 , KMAX - 1 
DO 778 I = 2 , [МАХ - 1 
778 DELTAT(I,K) = 1.0 
END IF 
ENO IF 


Cese THE DISSIPATION TERMS ARE CONSTRUCTED AND STORED IN THE ARRAYS DQ1, 
Csee 002,003 АМО 004. 
C 
CALL DISSIP 
C 
Cese THE RIGHT HAND SIDE AT KNOWN TIME LEVEL IS NOW COMPUTED AND ADDED 
CALL RESI (OMEGA) 
C 
Севе IF VISCOUS FLOW IS COMPUTED THE VISCOUS TERMS ARE ADDED TO DQ1 ETC. HERE. 
С 
IF(REYREF.GT.0.) CALL STRESS( ITN,DALFA) 
Се  [-SWEEP. 
C 
DTH = DT e 0.5 
DTW = DT « WWI 


DO 3 K - 2 , KMI 

CALL AMAT1(K,IMAX-1, XIX, XIZ, XIT) 

DO 4 I1 =1, 4 

00 4 12 = 1,4 

DO 5 1 = 2 , IMAX — 1 

Ее етк) = рш) ¢ DTH « DELTAT(I,K 

00(11,12,1-1,К) = — A(I1,12, 1-1) e OTH ° RA 
5 CONTINUE 
4 CONTINUE 
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Cees 


С 


ооо 


ООО 


O O O 


999 


991 


16 


17 
13 


IMPLICIT DAMPING ADDED HERE. 


DO 6 I1 ш1 , 4 

DO 6 I = 2 , IMAX - 1 

DT1 - DTW / YACOB(I,K) * OELTAT(I,K) 

Ge ME - Mm - 011 » YACOB( I-1,K) 
ЕЕ(11,11,1-1,К) « ЕЕ(11,11,1-1,К) - OT1 * YACOB(I*1,K) 
ММ(11,11,1-1,К) а 1. + 2. * DTW e DELTAT(I,K) 
CONTINUE 

DO 990 I - 2 , [MAX - 1 

GG(1,I-1,K) = DQ1(I,K) e DELTAT(I.K) 

2 та) - БОЗ К, » ОЕСТАТ p 

GG(3,I-1,K = 0Q3(1,K) * DELTAT(I,K 

GG(4,1-1,K) = 004(1,К) « ОЕЦТАТ(І,К) 

CONTINUE 

CONT INUE 


PERFORM BLOCK TRIDIAGONAL MATRIX INVERSION FOR THE ENTIRE PLANE 


CALL MATRX1( IMAX,KMAX) 


DO 991 K - 2 , KMAX - 1 
DO 991 I = 2 , ІМІ 
DQ1(I.K ш 00(1,1-1,К 
DQ2( I.K = 00(2,1-1,К 
DQ3(I,K ш 00(3,1-1,К 
DQ4(I.K = GG(4,I-1,K) 
CONTINUE 


K-SWEEP BEGINS HERE. 


DO 13 I -2 , IMI 

CALL AMAT2(I,KMAX-1,ZETAX,ZETAZ , ZETAT) 

00 15 11 =1, 4 

DO 15 12 301 ,54 

00 15 K = 2 , КМАХ - 1 

Бега ee - A(I1,I2, K*1)eDTH * OELTAT(I,K) 
00(11,12,1,К-1) = -A(11,12,K-1)*+DTH e DELTAT(1,K) 
CONT INUE 


SECOND ORDER DAMPING ADDED HERE. 


DO 16 I1 =2 1,4 

DO 16 K = 2 , КМАХ - 1 

DT1 = DTW / YACOB(I,K) * DELTAT(I,K) 
00(11,11,1,К-1) = DD(11,11,1,K-1) — DT! • Pron ER 
ЕЕ(11,11,1,К-1) ж ЕЕ(І1,11,1,К-1) - ОТ! + YACOB(1,K+1 
MM(11,11,1,K-1) = 1. + 2. e DTW e DELTAT(1,K) 

DO 17 K - 2 , КМАХ - 1 

GG(1,1,K-1 = 0Q1(I,K) 

00(2,1,К-1 - ai 

GG(3,1,K-1 = DQI(I,K 

GG(4,1,K-1 = DQ4(1,K) 

CONTINUE 

CONTINUE 


PERFORM BLOCK TRIDIAGONAL MATRIX INVERSION FOR THE ENTIRE PLANE 


CALL MATRX2( IMAX , KMAX) 
DO 18 I = 2 , IMAX — 1 
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DO 18 K = 2 , KM1 
DQ1(1,K) = GG(1,1,K-1) 
DQ2 a = GG(2,1,K-1) 
DQ3(1,K = GG(3,1,K-1) 
DQ4(1,K) з 00(4,1,К-1) 
18 CONTINUE 
E 
Cese UPDATE FLOW VARIABLES AT INTERIOR POINTS. 
967 CONTINUE 
RMAX = 0 
RUMAX = 0. 
БУМАХ ч 0. 
ЕМАХ = 0 
00 995 К ш 2 , КМ 
00 19 1 -» 2 , IMI 
Q1(I,K) = Q1(1,K) * DQ1(I,K) e сл ph 
к, - Q2(I,K) * DQ2(I,K) + YACOB(1,K) 
Q3(1,K) = Q3(1,K) + DQ3(I,K) * YACOB(I,K) 
Q4(1,K) ж 04(І,К) + DQ4(1,K) * YACOB(I,K) 
19 CONTINUE 
Cl 1! ! IDETERMINE WHERE IN FLOW FIELD DENSITY 1S CHANGING MOST RAPIDLY 
DO 995 1 = 2 , IMAX - 1 
IF (RMAX.LT. ass (001 (1. K)*YACOB(I,K))) THEN 
IR 
KR = k 
END IF 
RMAX z AMAX1(RMAX,ABS(DQ1(I,K) * YACOB(I,K)) 
RUMAX = AMAX1(RUMAX, АВ5(002(1. K) œ асов, d S 
RVMAX - AMAX1(RVMAX,ABS(DQ3(I,K) ev YACOB(I,K 
EMAX z AMAX1(EMAX,ABS(DQ4(I, K) e YACOB( I,K))) 
995 CONTINUE 
С ІЕ((ІТМ-1)/100%100.Е0.(ІТМ-1)) МКІТЕ (6,3002) 
IF(ITN .EQ. Ө) WRITE (6,3002) 
CIIIIISELECT INTERVAL AT WHICH OUTPUT OF RESIDUALS IS DESIRED 


O O O O 


IF((ITN-1)/10«10. EQ. ( ITN—-1)) WRITE (6,3001) RMAX,RUMAX,RVMAX , 
1EMAX , IR, KR 
RETURN 


3002 FORMAT(//,4X, 'DRMAX', 11X, 'DUMAX' , 11X, 'DVMAX' , 11X, 'DEMAX' , 10X, 


1' IR' , 3X, ' KR") 


3001 FORMAT(4(E14.8,2X),215) 


END 

SUBROUTINE MATRX 1 ( IMAX , KMAX ) 
COMMON/TRID/DD(4,4,161,41),MM(4,4,161,41), EE(4,4,161,41), 
166(4,161,41) 

COMMON/SCRAT/A(4,4,161) ,HH(4,4,161,41),C(4,5,161) 

REAL MM 

REAL L11,L21,L31,L41,L22,L32,L42,L33, LL 43, L44 


2 UFI, L21. C31, L41 


THIS SUBROUTINE PERFORMS THE BLOCK TRIDIAGONAL MATRIX IVERSION FOR 
AN ENTIRE PLANE DURING THE XI- SWEEP 


ОО 1 11 «1, 4 
DO 1 K = 2 , KMAX - 1 
АТ = 1. / MM(1,1,1,K) 


mm = GG(11,1,K) * А] 

НН(11,1,1,К) ж ЕЕ(11,1,1,К) * AI 
HH HIER - ра e Al 
ees Vek = EE(11,3,1,K) e Al 
HH(11,4,1,K) = EE(11,4,1,K) © AI 
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1 CONTINUE 


DD(I1,1, D; K) * GG 21 
00(11,2,1,К) ® 66(2,1-1,К 
DD(11,3,1,K) « GG(3,1-1,K) 
00(11,4,1,К) e GG(4,I-1,K) 


= ы e HH(1,12,1-1,K) 
- DO(11,2,1,K) * HH MU 
- DD rix) e НН(3,12,1-1,К 
- 00(11,4,1,к) « НН(4,12,1-1.К) 


3) © 21 
4) «121 


= L32 з U23 


4 — 132 • 024) » 151 
- 142 • 025 
(42 » 024 - (43 » U34 


C(1,1,K)) * L2I 

C(1,1,K) 

С(2,1,к)) < ШГ 

С(1,1,К) — 142 «СО ШЫ 
- (43 * С(3,1,К)) * L4I 


ELIO) + Ж 

G0. 2 h) 

С(2,2,К)) «131 

C(1,2,K). —- САКОО К) 
— L43 + C(3,2,K)) * L4I 


e C(1,3,K)) е 121 

* C(1,3,K 

* C(2,3,K)) * 31 

* C(1,3,K) — L42 * C(2,3,K) 
= 143» С(3,5,К)) е 141 


* C(1,4,K)) + L2I 


DO 1000 | = 2 , IMAX - 2 
DO 5 11441 , 4 
ОО 2 K -2 , KMAX - 1 
C(11,1,K) = GG(I1,I,K) — 
1 — 
2 ЯВ 
3 == 
2 CONTINUE 
DO 5 [2 м1, € 
DOS К =2, KMAX - 1 
А(11,12,К)= ММ(11,12,1,К) 
1 
2 
3 
С(11,1241,К)»- ЕНр 111211) 
5 CONTINUE 
DO 3 K = 2 , KMAX - 1 
L11 = A(1,1,K) 
СТ Гоа ТИЕ] 
012 = A(1,2,K) * LII 
U13 = A(1,3,K) * L1I 
014 = А Шал) e L1I 
L21 = A(2,1,K 
(31 = A(3,1,K) 
(41 « А(4,1,К) 
(22 = А(2,2,К) - 121 • 012 
L2I = 1. / L22 
U23 = (A(2,3,K) - L21 «e Un 
U24 = (A(2,4,K) — L21 + Un 
L32 = А(3,2,К) - 131 » U12 
(42 = А(4,2,К) – 141 • 012 
133 » A(3,3,K) - L31 • U13 
L3I = 1. / L33 
034 = (A(3,4,K) — L31 e U1 
L43 = ак - (41 » 1015 
(44 « А(4,4,К) - 141 • 014 
(41 "1. / 144 
21) = Clit Kk) e Е 
С(2,1,К) «ч (C(2,1,K) - 121 
C(3,1,K) = (C(3,1,K) — L31 
1 - L32 
C(4,1,K) « (С(4,1,К) - 141 
1 
со) ш ((1,2:К) «wu 
C(2,2,K) = (C(2,2,K) — L21 
C( (3.2 K) = (C (3.2 K) =S?Ë3] 
1 - L32 
C(4,2,K) = (C(4,2,K) - L41 
1 
С(1,3,К) = С(1,3,К) «111 
220) - зєв - 121 
С(3,3,К) = (С(3,3,Ку--131 
1 - 132 
С(4,3,К) « (С(4,3,К) - 141 
1 
С(1,4,К) * C(1,4,K) «111 
С(2,4,К) = СЕЛТ - 121 
С(3,4,К) «ч (С(3,4,К) - 131 


С(1,4,К) 
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ооо 


1 51132 «(24 Ку) 4151 
C(4,4,K) - (С(4,4,К) - 141 « С(1,4,К) - 142 « С(2,4,К) 
1 - 143 « С(3,4,К)) « 141 
С = C(1,5,K) * Lil 
C(2,5,K) = (C(2,5,K) - L21 * C(1,5,K)) * L2I 
C(3,5,K) » (C(3,5,K) - L31 * C(1,5.K) 
1 - 152 « ею . 131 
С(4,5,К) - (С(4,5,К) - 141 « С(1,5,К) - 142 « С(2,5,К) 
1 - 145 « С(3,5,К)) « 141 
С(3,1,К) = CRM - 034 « C(4,1,K) 
C(2,1,K) я С(2,1,К) - 924 » С(4,1,К) 
1 - 023 « de 1,K) 
C(1,1,K) «4 C(1,1,K) - 014 « С(4,1,К 
1 - 113» 1: 1 x - U12 * C(2,1,K) 
220) - Ээ) - 134 e 2,К 
C(2,2,K) = C(2,2,K) — U24 « 64.2.4) 
1 - 023 * C(3,2,K) 
C(1,2,K) » C(1,2,K) - 014 « С(4,2,К) 
1 - U13 * C(3,2,K) - U12 * C(2,2,K) 
C(3,3,K) 2 C(3,3,K) - U34 * C(4,3,K) 
C(2,3,K) » C(2, 3,K) - U24 * C(4,3,K) 
1 - 023 e C(3,3,K) 
C(1,3,K) - С(1,3,К) - Ul4 e C(4,3,K) 
1 - 013 « ак, - U12 * C(2,3,K) 
С(3,4,К) = С(3,4,К) - 034 » С(4,4,К 
С(2,4,К) = С(2,4,К) - U24 * C(4,4,K) 
1 — 423 « C(3,4,K) 
С(1,4,К) = С(1,4,К) - 014 * C(4,4,K) 
1 - 013. С(3,4,К) - 012 «е С(2,4,К) 
С(3,5,К) - С(3,5,К) - 034 « С(4,5,К) 
C(2,5,K) = C(2,5,K) - U24 * C(4,5,K) 
1 - U23 * C(3,5,K) 
C(1,5,K) = С(1,5,К) - 114 » С(4,5,К) 
1 — 413 * C(3,5,K) - U12 * C(2,5,K) 
3 CONTINUE 
DO 6 I1 21,4 
DO 9 K = 2 , KMAX - 1 
9 GG(I1,I,K) = C(I1,1,K) 
DO 6 I2 =1, 4 
DO 6 K = 2 , KMAX - 1 
НН(11,12,1,К) = С(11,12-41,К) 
6 CONTINUE 
1000 CONTINUE 
BACKWARD SUBSTITUTION 
00 7 1 = [MAX - 3, 1 , - 1 
DO 7 11 = 1,4 
DO 7 K - 2 , KMAX - 1 
66(11,1,К) » GG(I1,I,K) - HM(I1,1,I,K) e е 
1 - НН(11,2,1,К) » 00(2,141,К 
2 - НН(11.3,1,К) e GG(3,1+1,K) 
3 — HH(I1,4,I,K) e 66(4,1+1,К) 
7 CONTINUE 
RETURN 
END 
SUBROUT INE MATRX2( IMAX , KMAX ) 
COMMON/TRID/0D(4,4,161,41),MM(4,4,161,41), EE(4,4,161,41), 


166(4, 161,41) 
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COMMON/SCRAT/A(4,4,161),HH(4,4,161,41),C(4,5, 161) 
REAL MM 

КЕЛІ. (11,121,127 641 ,122,132 1421352 1432 14 
2,111,12Г 5 51:24) 


THIS SUBROUTINE PERFORMS THE BLOCK TRIOIAGONAL MATRIX IVERSION FOR 
AN ENTIRE J=CONSTANT PLANE DURING THE ZETA- SWEEP 


DO 1 I1 = 1, 4 
0011 = 2 , IMAX — 1 
AI ш 1. / ММ(1,1,1,1) 
66(11,1,1) = 66(11,1,1) * AI 
HH(11,1,1,1) = EE(11,1,1,1) * AI 
НН(11,2,1,1) = EE(11,2,1,1) © Al 
НН(11,3,1,1) = EE(11,3,1,1) © Al 
НН(11,4,1,1) = EE(11,4,1,1) * AI 
1 CONTINUE 
DO 1009 K = 2 , KMAX — 2 
DO 5 11= 1, 4 
00 2 | «2, 1МАХ- 1 
С(11,1,1) = 00(11,1,К) - 00(11,1.,1,К) ө 00(1,1,К-1 
1 - 00(11,2.,1,К) « 06(2.,1,К-1 
2 - орон! ° к 
3 — DD(I1,4,I,K) * GG(4,I.K-1 
2 CONTINUE 
DOS 12=1, 4 
005 І 22, IMAX - 1 
A(11,12,1)= MM(11,12,1,K) —- 00(11,1,1,К) * HH(1,I2, I, K-1) 
1 - 00(11,2,1,К) » НН(2,12,1,К-1) 
2 - 00(11,3.,1,к) * HH и 
3 - 00(11,4,1,К) « НН(4,12,1,К-1 
С(11,12-1,1)ж ЁЕ(11,12,1,Ю) 
5 CONTINUE 
DO 5122, IMAX - 1 
(11 ы А(1,1,1) 
Lilet. / L11 
U12 » A(1,2,I) * L1I 
013 = а e Lil 
014 = A(1,4,1) 4 111 
L21 = A(2,1.I) 
L31 us -4(3. 1,1) 
(41 = MA 
(22 т А(2,2,1) — L21 * U12 
L2Id» 1. / (22 
U23 = E - 121. 013) • 121 
024 « (А(2,4,1) - 121» па e L2I 
L32 = A(3,2,1) — L31 12 
(42 « А(4,2,1) - 141 « U12 
133 « А(3,3,1) - 151» 013 - 132 • U23 
L3I = 1. / L33 
0354 « (А(3,4,1) - 151 • 014 – 132 • 024) * L3I 
(43 ш ХОЛ - 141» 013 - 142» U23 
(44 = А(4,4,1) – 141 • 014 – 142 • U24 – 143 • U34 
141 ж 1, / L44 
SE = C(1,1,I) seg 
С(2,1,1) т (0С(2.1,1) -121 «С(1:1:1)) 220221 
C(3,1,1). «e (C(3,1.1) - 2 ЇМ 
1 - (52. CD “ |51 
С(4,1,1) ы» (С(4,1,1) - 141 « С(1,1,1) - 142» С(2,1,1) 
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GE Со уе 
С(2,2.1) = (С(2,2,1) — 121 e 
C(3,2,I) = (C(3,2,I) - (31 9 
— L32 * 
С(4,2,1) « (С(4,2,1) - 141 * 
C(1,3,I) = C(1,3,I) * Lil 
C(2,5 1) ча (2,3,17”-- 12154 
C(3. 3.1) = (C(3,3.1) =- L31 = 
- 132 ® 
С(4,3,1) » (С(4,3,1) - 141 » 
C(1,4,I) = C(1,4,1) * L1I 
С(2,4,1) = (С(2,4,1) - 121» 
C(3,4,1) *»»(C(5,4,I) - 151 е 
= 132 • 
C(4,4,1) = (C(4,4,1) – 141 • 
СО 5.1) зот, ГЕ! 
С(2.5.1) аВ(С(2.5.,1) - 121 , 
c(3.5.1) ¿mr (0(3,5,1) — L31 = 
= 132 » 
C(4,5,1) = (C(4,5,1) - 141 » 
1 
С(3,1,1) = С(3,1,Г) - 034» 
С(2,1,1) « С(2,1,!) - 424 * 
1 - (23 » 
C(1,1,I) = C(1,1,1) - 014% 
- (15 e 
c = С(3,2,1) - U34 « 
С(2,2,1) ж С(2,2,1) - 024» 
- 23 » 
С(1,2,1) « С(1,2,1) - 014 • 
= 13 » 
C(3,3,1) = C(3,3,1) - 034 » 
С(2,3,1) = C(2,3,1) - U24 * 
= 023 г, 
ССТ З ID) 2wC(1,3,I) —- Uf4 e 
= 13 » 
С(3,4,Г) = С(3,4,Г) - 054» 
С(2,4,1) = C(2,4,1) - U24 * 
= 023 • 
C(1,4,I) = С(1, 4,1) – 014 » 
1 = 013 • 
С(3,5,1) = C(3,5,1) - U34 • 
С(2,5,1) м С(2,5,1) - 024 * 
- 23 » 
С(1,5,1) » C(1,5,I) – 014 s 
- U13 s 
CONTINUE 
DO 6 I1 = |, 4 
DO 9 I = 2 , [MAX - 1 
GG(I1,I,K) = ((11,1,1) 
DO 6 I2 = 1 , 4 
DO 6 I = 2 , [MAX - 1 
HH(11,12,1,K) = C(11,12+1,1) 
CONTINUE 
CONT INUE 


2394374 C(3,1,[)) * L4I 


002227 

С(1,2.1) - 142 « С(2,2,1) 
- 145 » С(5,2,1)) « 141 

AE . 121 

СОТ) 

2(2:321)) 4131 


C(1,3,1) — L42 < C(2,3,I) 
-143 * C(3,3,1)) е“ 141 


— L43 « С(3,4,1)) » 141 


С(1:522)) 742121 

ШІ) 

C(2,5,1)) * 131 

С(1,5,1) - 142» С(2,5,1) 

- 143 • С(3,5,1)) 4 141 

м 
с(4,1,1) 
C(3,1,I) 
С(4,1,1) 
C(3 3.1) 
С(4,2,1) 
C(4,2,1) 


G2 


U12 * 


62.2, 1) 


- 012 » С(2,5,1) 


O 
~ 
CA 


012 » С(2,4,1) 


€(2,5,1) 
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BACKWARD SUBSTITUTION 


ООО 


DO 7 K = KMAX - 3, 1 , - 1 

DO 7 11 = 1, 4 

DO 7 I = 2 , [МАХ - 1 

GG(11,1,K) = GG(11,1,K) - HH(11,1,1,K) e GG(1,1,K+1) 
1 - HH(11,2,1,K) * GG Son 
2 - HH(11,3,1,K) * GG(3, I, K*1 
3 - HH(11,4,1,K) e GG(4,1,K+1) 
7 CONTINUE 

RETURN 

END 

SUBROUT INE METRIC 

COMMON/F IX/OMEGA 


COMMON/DGRID/OT , IMAX,KMAX, I TEL, I TEU 
COMMON/GRID1/X(161,41),2(161,41) 

COMMON/GRID/YACOB(161,41) 
COMMON/MTRIX/XIX(161,41),XIZ(161,41), ZETAX(161,41), ZETAZ(161,41), 
1XIT(161,41),ZETAT(161,41) 


C 

Coro SUBROUTINE METRIC COMPUTES THE METRICS IN BOTH DIRECTIONS AND 
C THE UNSTEADY COEFFICIENTS ETAT, ETC. 

C 


DO 1000 K = 1 , KMAX 
DO 1000 I = 1 , IMAX 
XTAU = OMEGA e Z(I,K) 
YTAU = OMEGA + (-X(1,K)) 
Севе PRESENT SET UP IS FOR FLOW PAST AN AIRFOIL. 


C 
CIIIIICENTRAL DIFFERENCES AT INTERIOR POINTS, TWO-POINT ONE-SIDED 
CI!!! !IDIFFERENCES DOWNSTREAM, THREE-POINT AT OTHER OUTER BOUNDARIES 
IF(I.EQ.1.0R.I.EQ. IMAX) GO TO 1e 
XXI з .5 * (X(I*1,K)-X(I-1,K)) 
ZXI = .5 *« (Z(I*1,K)-Z(I-1,K)) 
GO TO 15 
10 IF(I.EQ.IMAX) GO TO 16 
XXI = 1.9 e (X(2,K) - X(1,K)) 
2Х1 «1.0» (7(2,К) - Z(1,K)) 
СО ТО 15 
16 XXI = 1.9 e (X(IMAX,K) - X(IMAX-1,K)) 
2Х1 = 1.0 * (Z(IMAX,K) - Z(IMAX-1,K)) 
15 CONTINUE 
IF(K.EQ. 1. OR.K. EQ. KMAX) GO TO 17 
XZET = .5 о 
ZZET = .5 4(2(1,К-1)-2(1,К-1)) 
GO TO 20 
17 IF(K.EQ.KMAX) GO TO 18 
XZET*» 2. « X(I1,2)-1.5 « X(I,1) - .5 * X(I,3) 
ZZ "2, * 2(1[,2) =- 1.5. « ZCE, 1941555 2200055 
GO TO 20 
18 XZET = 1.5 e X(I,KMAX)-2.0+ EE 
22ЕТ = 1.5 ® 2(1,КМАХ)-2.9 2(1,КМАХ-1)4.592(1 ,КмАХ-2 
20 CONTINUE 
С!!! ! ! COMPUTE JACOBIAN 


YACOBI = XXI « ZZET - XZET e ZXI 
YACOB(I,K) - 1. / YACOBI 

XIX(I,K) = ZZET е ҮАСОВ(1,К) 
UNES = -XZET + YACOB(I,K) 
XTAU = OMEGA + Z(1,K) 
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YTAU = — OMEGA , Х(1,К) 

XIT(I,K) = — XIX(I,K) * XTAU - XIZ(I,K) * YTAU 

SE = —ZXI * YACOB(I,K) 

ZETAZ(1,K) = XXI + YACOB(I,K) 

ZETAT(I,K) » — ZETAX(I.K) * XTAU - ZETAZ(I,K) * YTAU 
CONT INUE 

RETURN 

END 

SUBROUTINE DISSIP 
COMMON/FLOW/Q1 (161,41) ,Q2(161,41) ,Q3(161,41) ,Q4(161, 41) 
COMMON/PERTR/0Q1(161,41) ,002(161,41) ,0Q3(161,41) ,0Q4(161,41) 
COMMON/DGRID/OT, IMAX,KMAX, ITEL, ITEU 
COMMON/GRID/YACOB(161,41) 

COMMON /DAMP /WW , WWI 

DIMENSION P(161),€PS(161),01S1(161,4) ,DIS2(161,4) 


THIS SUBROUTINE ADDS THE FOURTH ORDER DISSIPATION TERMS TO THE 
RIGHT HAND SIDE 


IMI = [МАХ — 1 
КМ1 = КМАХ - 1 
ІМ2 - ІМАХ - 2 
КМ2 = КМАХ - 2 


DO 10 K=2 , KM1 

COMPUTE SWTICHING FUNCTION BASED ON SECOND DERIVATIVE OF PRESSURE 
DO 1 I = 1 , IMAX 

P(I) 2 .4 * (Q4(I,K)-(Q2(I,K)**24Q3(I,K)**2)/(2. Q1 (I,K))) 
DO 2 I =1 , IMAX 

IP2 = I + 2 

IF(I.EQ.IM1) IP2 = IMAX 

IM2 = 1-2 

IF(I.EQ.2) IM2 = 1 

IP1 = I + 1 

IF(I.EQ. IMAX) IP1 2 [МАХ 

ЇМ = | — 1 

IF(I.EQ.1) IM = 1 

IF(I.EQ.1) IM2 = 1 

IF(I.EQ. IMAX) IP2 = IMAX 

ЕР5(1) » ABS(P(IP1)-2. *P( I) 4P( IM))/ABS(P(IP1)*2.*P( I) 4P( IM)) 
VOL = 2. /(YACOB(I,K)*YACOB( IP1,K)) 

VOL = 1. 

01$1(1,1) = (QI!(IP1,K)-Q1(I,K))*VOL 

0151(1,2) = (Q2(IP1,K)-Q2(I,K))*VOL 

0151(1,3) » (Q3(IP1,K)-Q3(I,K))*vOL 

HP1 = Q4(IP1,K)+P(IP1) 

HP = Q4(1,K)+P(1) 

НМ1 = 04(1М,К) 4 Р(1М) 

HP2 = Q4(1P2,K) + P(IP2) 

HPM * Q4( IM,K)*P( IM) 


01$1(1,4) = (HP1-HP)+*VOL 
01$2(1,1 Sede cU EU qud Ké -Q1(IM,K кү 
01$2(1,2 Q2(IP2,K)-3.*(Q2(IP1,K)-Q2(1,K))-Q2( IM.K)) «VOL 


m= 
a= 
a 
=== 


DIS2(I,3) (Q3(IP2,K)-3.«(Q3(IP1, K)-Q3( I, K) )-Q3( IM.K) ) «voL 
DIS2(1,4) (HP2-3.*(HP1-HP)-HPM) «VOL 
CONTINUE 


0015 141 , IMI 
D2P = AMAX1(EPS(I),EPS(I+1)) 
С22 = 60. • 02Р 

С11 = АМАХ1 (0.0, (1.-С22)) 


E 


C22 = C22 4 eh 1 
C11 = C11 + WW/YACOB(I,K) + DT 


CIITII I SWITCH ON SECOND-ORDER AND SWITCH OFF FOURTH-ORDER DISSIPATION 
CI!11! IN VICINITY OF SHOCKS 
01$1(1,1) = С11 • 0152(1,1) + С22 • 0151 n 
01$1(1,2) = С11 ® 01$2(1,2) + 622 * 01$1(1,2 
DIS1(1,3) = C11 • 0152(1,3) + C22 * DIS1(1,3) 
0151(1,4) = C11 e DIS2(1,4) + C22 « 0151(1,4) 
15 CONTINUE 
0016 І - 2, ІМІ 
001(1,К) « 0151(1,1) - 01$1(1-1,1 
002(Т.К) = 01$1(1,2) - 01$1(1-1,2 
003(Т,К) = 01$1(1,3) - 01$1(1-1,3 
004(Т,К) = 01$1(1,4) - 01$1(1-1,4 
16 CONTINUE 
10 CONTINUE 
C K OIRECTION 
C!!!t! FOURTH-ORDER DISSIPATION ONLY 
DO 30 I «2, ІМІ 
WT- 0.5 «€ DT « WW / YACOB(I,2) 
МЗ = 2.5 « ОТ « МИ / ҮАСОВ(1,КМ1) 
001(1,2) «ит» (01(1,1) - 2. « 01(1,2) + Q1(1,3)) 
1%001(1,2) 
002(1,2) зэМТ» (02(1,1) - 2. з 02(1,2) + 92(1,3)) 
1%002(1,2) 
003(1,2) -ЯТ« (03(1,1) - 2. « 03(1,2) + Q3(1,3)) 
1%003(1,2) 
004(1,2) «ЯТ» (04(1,1) - 2. » 04(1,2) + 04(1,3)) 
14004(1.2) 
WT= W3 
DQ1(I,KM1) =WT* (Q1(1,KM2) - 2. * Q1(I,KM1) + Q1(1,KMAX)) 
14001(1,КМ1) 
002(1,КМ1) ачит» (02(1,КМ2) - 2. » 02(1,КМ1) « 02(1,КМАХ)) 
1+DQ2(1,KM1) 
DO3(1,KM1) -WTs (Q3(I,KM2) - 2. © Q3(1,KM1) * Q3(I,KMAX)) 
14003(1,КМ1) 
DO4(I,KM1) -WTe (O4(I,KM2) - 2. * O4(I,KM1) tr Q4(I,KMAX)) 
1-004(1,КМ1) 
DO 35 K = 3 , KM2 
Wis — ОТ « WW / YACOB(I,K) 
DO1(1,K) =WT* (01(1,K+2) — 4. e Q1(1,K+1) + 6. © QI( 
11,K) — 4. «€ Q1(I,K-1) + 01(1,K-2))+D01(1,K) 
002(1,К) =WTe(Q2(1,K+2) - 4. © Q2(1,K+1) + 6. * Q2( 
11,К) - 4. «€ 02(1,K-1) + 02(1,K-2))+002(1,K) 
003(1,К) sWTe(Q3(I,K42) - 4. © Q3(1,K+1) + 6. * Q3( 
11,K) - 4. «€ 03(1,K-1) + Q3(1,K-2))+D03(1,K) 
004(1,К) «МТг(04(1,К-2) - 4. » 04(1,К41) + 6. • 04( 
11,К) - 4.» 094(1,К-1) + 04(1.К-2))4004(1,К) 
35 CONTINUE 
39 CONTINUE 
C 
RETURN 
ENO 


SUBROUT INE WALLBC 
COMMON/SURF /PSUR( 161) 
COMMON/GRID1/X(161,41),Z(161,41) 


COMMON/PAR/GAMMA ,REYREF ,ALFA,ALFA1 REDFRE,AMINF ,ALFAI 


COMMON/DGRID/DT , IMAX,KMAX, ITEL, ITEU 
COMMON/GRID/YACOB(161,41) 
COMMON /DAMP /WW , WW I 
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COMMON/MTRIX/XIX(161,41),XIZ(161,41) , ZETAX(161,41) , ZETAZ(161,41), 
1XIT(161,41), ZETAT(161, 41) 
COMMON/FLOW/Q1 (161,41) ,Q2(161, 41) ,Q3(161,41) ,Q4(161,41) 
DIMENSION C1(4) 
DIMENSION A(2,2),RHS(2) 
C!!!!!APPLY BOUNDARY CONDITIONS ON THE CUT AND THE AIRFOIL SURFACE 
DO 9 I=ITEU, IMAX 
11 = [МАХ + 1 -I 
Q1(I,1) = .5 +. (Q1(1,2)+01(11,2)) 
02(1,1) =.5•(02(1,2)+02(11,2)) 
03(1,1) = .5 • (81122) % 03(11,2)) 
04(1,1) « .5 » (04(1,2)404(11,2)) 
01(11,1 SM 
Q2 e 1,1 
Q3(11,1)2Q3(I,1) 
Q4(I11,1)2Q4(I, 1) 
9 CONTINUE 
00 1 I= ITEL , ITEU 
К = 3 
ea -» XIT(I.K) 
C1 J - цо 
C1(3) = XIZ(I,K 
UCON3 » (Q2(I,K)*C1(2) 4Q3(I,K)*C1(3)) 
1/Q1(1,K) 
K = 2 
C1(1) 
d 
C1(3 
UCON2 
1⁄Q1(I,K) 
RHS(1) = 2. * UCON2 - UCON3 - XIT(I,1) 
C FOR VISCOUS FLOWS SET UCON TO ZERO ALSO 
IF(REYREF.GT.0.) RHS(1) = - XIT(I,1) 
А(1,1) = сш 
A(1,2) = XIZ(1,1) 
А(2,1) = 2ЕТАХ(І, 1) 
A(2,2) = ZETAZ(I,1) 
RHS(2) = - ZETAT(1,1) 
TEMP1 = Ke 
TEMP2 = А(1,2) 
ТЕМРЗ = 220) 
ТЕМР4 = А(2,2) 
DEN = 1. /(TEMP1 e TEMP4 — TEMP2 • ТЕМРЗ) 
A(1,1) = A(2,2) * DEN 
A(1,2) 2 - TEMP2 • DEN 
Au - — TEMP3 • ОЕМ 
A(2,2) = TEMP1 * DEN 
SAU. = 2. 4 01(1,2) - 01(1,3) 
Q2(1,1) e Q1(I,1)*(A(1, 1) *RHS(1)*A(1,2) *RHS(2) 
Q3(1,1) = Q1(1,1) о 
CONT INUE 
DO 19 I=mITEL ,ITEU 
они 
W2=Q3(1,2)/01(1,2) 
Р25(САММА-1.)4(04(1,2)-9.5401(1,2)4(026024М242)) 
03-02(1,3)/01(1,3) 
W32Q3(1,3)/Q1(1,3) 
P3=(GAMMA—1.)+(04(1,3)-0.5:Q01(1,3)+(U3eU3+W3+W3)) 
Р1-(4.»Р2-Р3)/3. 
РУОВ (1 )=(САММА*«Р1-1.)/(.7*АМ1МЕ®*2) 


XIT(I.K) 

ХЇХ(1,К 

XIZ(I.K 
(Q2(1,K)*C1(2)4Q3(I,K)*C1(3)) 


—^ 
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19 


dl Тос 
W1=Q3(1,1)/01(1,1) 
04(1,1)4Р1/ (САММА-1. )40.5401(1,1)4(0164Л 90101) 

КЕТОВМ 

ЕМО 

SUBROUTINE STRESS(ITN,DALFA) 

COMMON/F LOW/Q1 (161,41),02(161,41),03(161,41),04(161,41) 
COMMON/OGRID/OT , IMAX, KMAX, ITEL, ITEU 
COMMON/GRID1/X(161,41),2(161,41) 
COMMON/PAR/GAMMA , REYREF, ALFA, ALFAT, REDFRE,AMINF,ALFAI 

COMMON PERTE/DOI (161 =e) ) Doacne t spana M EDD анин 
COMMON /MUTUR/CMU( 161,41) 

DIMENSION AA(161, 41) 

1, RH1(161) ,RH2(161) , RH3(161) , RH4(161) 
COMMON/LOGIC/RSTRT,PITCH,PLUNGE 

LOGICAL RSTRT,PITCH,PLUNGE 

VOTED - Бос) / Q1(I,J) 

V(I1,J) = Q3(I,J) / Q1(1,J) 

THIS SUBROUTINE ADDS VISCOUS TERMS TO THE RIGHT HAND SIDE 
GOGM = GAWA / (GAMMA — 1.) 

IF(ITN.GT.18.OR.(RSTRT)) CALL EDDY(DALFA) 

COMPUTE U AND 
KMAXM1 = КМАХ 
ЇМАХМ1 = [МАХ 
PR = 1. 

00 19 К 21 , KMAX 

00 10 I 2 1 , IMAX 

E = Q4(1,K) / QI(I,K) - 0.5 * (U(I,K)**2*V(I,K)*9*2) 


1 
1 


I I < 


10 AA(I,K) « DON a E 


COMPUTE TXX,TXY AND VISCOUS DISSIPATION AT I = 1 / 2 


DO 30 K = 2 , KMAXM1 

DO 28 I = 2 , IMAX 

UXI 2 U(I.K) - U(I-1,K) 

VXI = V(I,K) - У(1-1,к) 

АХ] = AA(I,K) — АА(1-1,К) 

UZET= E 
VZETe .25«(V(I,K41)-V( I, K-1)«V(I-1,K*1)-V( I-1,K-1 
AZET= .25*(AA(I,K+1)-AA(I,K-1)+AA(I-1,K+1)-AA(I-1,K-1)) 
ХХ| - 2 - X(I-1,K) 

ZXI = Z(I,K) -2(1-1,К) 

XZET= .25 e QI KD XC Kendo Л е Т) 
ZZETe .25 , (2(1,К41)-2(1,К-1)-2(1-1,К41)-2(1-1,К-1) 
УАС = ХХ! ® ZZET - ZXI e XZET 

YAC = 1. / YAC 

XIX = ZZET e YAC 

2ЕТАХе - 2Х| « ҮАС 

XIZ = -XZET e YAC 

ZETAZ= XXI * YAC 

CNM = .5 * (CMU(I,K) + CMU(I-1,K)) 


UX = UXI $ XIX + UZET » ZETAX 
VX = VXI * XIX 4 VZET e ZETAX 
АХ = AXI * XIX + А2ЕТ * ZETAX 
UZ = UXI * XIZ + UZET * ZETAZ 
М2 = VXI e XIZ + VZET * ZETAZ 
AZ = AXI * XIZ + AZET * ZETAZ 


ТХХ = -(-4. * UX + 2. * VZ) ,« ОМ / 3. 
TXY = CNM * (UZ + VX) 
ТҮҮ = —CNM / 3. e (-4. < VZ + 2. © UX) 
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54 = ((U(I,K)+U 
1 


+ CNM / PR/(GAMMA – 1.) • 
1-1 х: í. -1,K))+TYY)+0.5 
+ CNM / PR / (GAMMA — 1.) © AZ 


R4 » ((U(I, сл -1 LO CURA. K)+V(I-1,K))+TXY)+0.5 


DEBUG 

TURN OFF ENRGY DISSIPATION AND DIFFUSION 
R4 = 0. 
S4 = 0. 
о, = 0. 
RH2(1) = (XIX • ТХХ + XIZ » TXY) / YAC 
а, = ыг • TXY + ХІ2 • ТҮҮ) / ТАС 
БН4(1) = (ХІХ * R4 + XIZ * S4) / ТАС 


DO 30 1 «2, ІМАХМ1 


001(1,К) «ч 001(1,К) + RH1(I+1) - RH1(1) 
Si? = DQ2(1,K) + RH2(I+1) - atn 
003(1,К) « 003(1,К) 4 RH3(I41) - RH3(I) 
DQ4(I,K) = DQ4(I,K) + RH4(I+1) — RH4(I) 
CONTINUE 


IN THE Z DIRECTION 
DO 70 [ = 2 , ІМАХМЇ 


DO 60 K = 2 , KMAX 
UXI = .25 * (U(I+1,K)-U(I—-1,K)+U(I+1,K—- Dus ШЕ 1,К-1) )) 

VXI = 28  (У(141,К)-Ч(1-1,К)4-У(141,К-1)-4(1-1,К-1) 

AXI = .25 « (АА(141,К)-АА(1-1,К)4АА(141,К-1)-АА(1-1,К-1)) 
XXI = .25 » И Dna -1, apex. ге ЇЕ 1,К- H 

ZXI = .25 « (Z(I41,K)-2(I-1,K)42(I41,K-1 1-1,К-1) 


UZET = ЗИ - U(I,K-1) 
У2ЕТ « У(1,К) - У(1,К-1) 

А2ЕТ - АА(1,К) - АА(1,К-1) 
XZET = X(I,K) —- X(I,K-1) 

ZZET = Z(1,K) - Z(1,k-1) 

ҮАС = ХХІ • 27ЕТ — ZXI * XZET 
ҮАС = 1. / ХАС 

XIX = ZZET * УАС 

ZETAXe - ZXI e YAC 

Х12 « -Х2ЕТ + YAC 

ZETAZ= XXI * YAC 


CNM = ‚5 * (CMU(I,K) + CMU(I,K-1)) 

UX = UXI * XIX + UZET * ZETAX 

VX = VXI * XIX + VZET • 2ЕТАХ 

AX = AXI * XIX + AZET * ZETAX 

UZ = UXI * XIZ + UZET * ZETAZ 

VZ = VXI * XIZ + VZET * ZETAZ 

А2 = AXI * XIZ + AZET * ZETAZ 

TXX = —(-4, © UX + 2. * VZ) * CMM / 3. 
TXY = CNM • (02 + VX) 

TYY = -CNM / 3. e (-4. e VZ + 2. as UX) 
R4 » ((U(I,K)*U(I,K-1)) eTXX*(V(I, K)*V(I,K-1)) 9 TXY)«0.5 


4 CNM / PR/(GAMMA — 1.) e AX 
54 - ((U(I,K)*U(I,K-1)) TXY-(V(I, Qua. K-1))+TYY)*0.5 
+ СММ / PR / (САМА — 1.) » 


R4 = 0. 

S4 = 0. 

RH1(K) = 0. 

ЕН2(К) - (?ЕТАХ » ТХХ + ZETAZ + TXY) / YAC 
RH3 к = Е * TXY + ZETAZ «€ TYY) / YAC 
RH4(K) = (ZETAX * R4 + ZETAZ » S4) / YAC 


DO 70 K 2 2 , КМАХМ1 
001(1,К) » DQ1(I,K) 4 RH1(K41) - RH1(K) 
002(1,К) « 002(1,К) 4 RH2(K41) - RH2(K) 
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70 


400 


БТ = DQI(I,K) + но - и, 
094(1,К) = 004(1,К) + RH4(K+1) - ВН4(К 
CONT INUE 


RETURN 

END 

SUBROUTINE LOAD(CPS,CL,CD,CM,ALFAS) 
COMMON/GRID1/X(161,41),Y(161, 41) 
COMMON/DGRID/OT, IMAX,KMAX, ITEL, ITEU 
DIMENSION CPS(161) 


THIS SUBROUTINE COMPUTES THE INVISCIO CONTRIBUTIONS 
TO LOADS ON THE AIRFOIL SURFACE 


IMAXM1 з= [МАХ — 1 


CL = @. 

СО = Q. 

CM = 0. 

00 400 І = ITEL , ITEU - 1 
XL = .5 « ИЕ 
УЕ а .5 4 (Ү(Т,1)4Ү(141,1)) 
ОХ = Х(1+1,1) - X(I,1) 


OY = Ү(141,1) - Ү(1,1) 

CPA = CPS(I+1) « .5 + CPS(I) « .5 
OCL = CPA e (-DX) 

OCD = CPA « DY 

CL = CL + DCL 

CO = CD + DCD 

CM = CM + DCD • YL — OCL + XL 


DCL = CL + COS(ALFAS) — CD e SIN(ALFAS) 

CD = CL + SIN(ALFAS) + СО • COS(ALFAS) 

CL = DCL 

RETURN 

END 

SUBROUTINE WRAP(II,JJ,XSING, YSING, XP, YP, SO, A0 , YO) 
DIMENSION S0(161,4),Y0(41,4),A0(161,4) ,XP( 1) . YP(1) 


THIS SUBROUTINE UNWRAPS THE AIRFOIL AND STORES THE UNWRAPPED 
SURFACE IN ARRAYS A@ AND S@. IT ALSO DETERMINES THE STRETCHING 
IN THE ETA DIRECTION. 


IMIO = (II + 1) / 2 

DY = .8 / (JJ - 2) 

001 J =2 , JJ 

Y = FLOAT(J-2) * DY 

YO(J,1) = 1.25.» Y / (1. — Y ° Y) 
YO0(1 , 1) = — Ү9(3,1) 

PI = 4. • ATAN ( 1.) 

ANGL = PI + PI 

U = XP(1) - XSING 

V = YP(1) - YSING 

Ц » 1. 

V = 0. 

IIM1 = II - 1 

DO 2 I - 1 , II 

X11 = XP(I) - XSING 

Ү11 = YP(I) - YSING 

ANGL = ANGL + ATAN2((UsY11-VeX11),(UsX11+VeY11)) 
R = SQRT(X11ee2 + Y11e 2) 


U e Vi? 
V =Y11 

R = SQRT(R) 

Ад(1,1) 4 8 » COS(.5 e ANGL) 


2 SO(1,1) = R + SIN(.5 e ANGL) 
CI!!1! IF OUTPUT OF UNNRAPPED COORDINATES IS DESIRED 
C WRITE (5,1000) 
С WRITE (6,2000) (1,A0(1,1),S0(1,1),1 = 1 , 11) 
RETURN 
1000 E тогвЕр COORDINATES IN THE TRANSFORMED PLANE’ ) 
2000 FORMAT(I5 , 2F20.8) 
END 
SUBROUTINE TABINT(XP,YP,XSING,YSING,N) 
DIMENSION XP(161), YP(161),S0(161),A0(161) 
C!! II !SMOOTH THE AIRFOIL SURFACE BY FINDING ADDITIONAL POINTS 
U - XP(1) - XSING 
V = YP(1) - YSING 
U x 1. 
У = 0. 
АМСЕ = 8. * АТАМ(1.) 
DO 1 I = 1,N 
X11 » XP(I) - XSING 
У11 = YP(I) — YSING 
ANGL = ANGL + ATAN2((UsY11-Vex11),(Usx11+VeY11)) 
R = SQRT(X11e9e2 + Y11 ее 2) 
U = X11 
V = Y11 
R = SQRT(R) 
AO(I) =R à у ‚ .5) 
1 50(1) = В * SIN(ANGL * .5) 
Dx =(AQ(N)-AQ(1))/96. 
АОЗТ = AQ(1) 
DO 3 I = 1 , 97 
XX = SNE ‚ ОХ + АОЗТ 
CALL TAINT(AO,SO,XX, YY, N, S, NER, MON) 
TOU = XX * XX — YY * YY + XSING 
3 YP(I) = 2. © XX s YY + YSING 
RETURN 
END 
SUBROUTINE TAINT(XTAB,FTAB,X,FX,N,K,NER,MON) 
DIMENSION XTAB(1),FTAB(1),T(10),C(10) 
C 
C NASA — AMES SUBROUTINE FOR POLYNOMIAL INTERPOLATION 
C OF A TABULATED FUNCT ION 
C 
IF(N-K) 1,1, 2 
1 NER = 2 
RETURN 
2 IF(K-9) 3,3,1 
3 оос 4,4,5 
5 ІҒ(МОМ-2) 6,7,4 
4 J =O 
NMÍ 2 N- 1 
DO 8 I = 1 =, ММ! 
IF(XTAB(I) - XTAB(I+1)) 9,11,10 
11 NER = 3 
RETURN 
9 Ј ж Ј-1 
GO TO 8 
10 J = 4+1 
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O O O O 


14 


13 


17 
16 
18 


19 
29 


22 
21 
26 


23 


25 


24 


CONT INUE 
MON = 1 
ІК(ТІ72 77767776 
МОМ ж 2 
DO 13 1 - 1 , N 
ЇЕ(Х - XTAB(I)) 14,14,13 
4 = [ 
GO TO 18 
CONTINUE 
GO TO 15 
DO 16 І з 1 , N 
IF(X-XTAB(I)) 16, 17, 17 
J = I 
GO TO 18 
CONTINUE 
J =N 
J = J — (км) / 2 
IF(J) 19,19,20 
4 = 1 
M= J + K 
IF(M - N) 21,21,22 
) =) — 1 
GO TO 20 
KPI = K + 
JSAVE = J 
DO 23 L = 1, KP1 
C(L) = X — XTAB(J) 
T(L) = FTAB(J) 
J = J+! 
DO 24 J = 1,K 
I = J+! 
E IEA 
= [4 
1Е(1-КР1) 25,25,24 
CONTINUE 
FX = T(KP1) 
NER = 1 
RETURN 
END 
SUBROUTINE SING(N2,N,X,Z,XLE,YLE,TEA,TES,XSING,YSING ,CHD) 


1 


THIS SUBROUTINE COMPUTES SINGULAR POINT LOCATIONS. 
DIMENSION X(2) . Z(2) 


NU = N2 

N1 = N2 + 1 

N3 2N2- 1 

X1 = AS 

71 2:7 13) 

Х2 = Х(М2 

22 = e. 

X3 -X NaS 

13 = Z(N3 

01 = X2 se 2 - XI e. 2 
02 = 22 ве 2 - 21 ef 2 
03 » X2- X1 

04 2 22 - 21 

05 = ХЗ ве 2 – ХЇ в• 2 
06 = 23 eo 2 - 21 •• 2 
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10 
20 


C 
2111! 
8 


1 


07 = Хз -Х! 
08 = 73 - 21 

B = үр: .(01-02)- 034(05406))/(2.4(07404-03408)) 
IF(ABS(D3).LT.ABS(D7)) GO TO 10 

А = (01 + D2 — 2. • В * D4) / (2. * 03) 

GO TO 2@ 

А = (05 + 06 – 2. • В • 08) / ( 2. * D7) 

CONT INUE 

В = SQRT((X2-A)»« 2 4 (22-8)ө2) 

XLE = X(NU) 

YLE = Z(NU) 

CHO = X(1) - X(NU) 

А2 = 2 2-10) /(X(2) - x(1)) 

А1 » (Z(N)-Z(N-1))/(X(N)-X(N-1)) 

TES = .5 ® (A1 4 A2) 

ТЕА = А2 - A1 

TEA = TEA * 57. 29578 

XSING = ЕЕ 12. 

YSING = (B+YLE) / 2. 

RETURN 

END 

SUBROUTINE AIRFOL(II,JJ, IT, IE) 
COMMON/GRID1/X(161,41),2(161,41) 

COMMON /YS YM/ I SYM 

DIMENSION S0e(161,4),A0(161,4),Y0(41,4), XP( 161) , YP( 161), 
1Е(161),Е(161),80(49) 


DATA (B0(1),121,32)/1.,1.0414,1.0836,1.1270,1. 1715, 1.2175, 1.2651, 
11.3145,1.3659,1.4199,1.4755,1.5349,1.5973,1.6636,1.7342,1.8099, 
21.8914,1.9799,2.0764,2.1829,2.3012,2.4341,2.5653,2.7597,2.9646, 
33.2106,3.5141,3.9019,4.4219,5.1687,6.3632,8.6809/ 


|СОМРИТЕ THE COMPUTATIONAL GRID POINTS BASED ON INPUT AIRFOIL SHAPE 
0081=1, 32 

Ад(1,1) = 80(1) 

READ (5,1) 

FORMAT ( 1X) 

READ (5,2) FNU,FNL,ZSYM 
FORMAT ( 3F 10.0) 

ISYM - 0 

IF(ZSYM.NE.@.) ISYM = 1 
11 = 157 

JJ =m 41 

IT = 31 


PI = 4. • АТАМ(1.) 

NU = FNU 

NL = FNL 

N = NU + NL - 1 

READ(5,1) 

READ (5,333) (XP(I), YP(1),I = NL , N) 
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C 


C 


333 
9994 


9995 


16 
9996 


9997 


6791 


438 
439 


ЕЕ 
FORMAT(F20.8) 

L=N +1 

IF(ZSYM .GT. Ө.) GO TO 9995 

L = NL + 1 

READ(5, 1) 

READ (5,333) (XP(L-I), YP(L-I) , I1, NL) 
GO TO 9996 

K1 = L 

DO 16 I NL, N 

K = Кі — Í 

XP(K) = XP(I) 

УР(К) = - ҮР(1) 

CONT INUE 

SCALE = 1. /(XP(1)-XP(NL)) 

XX = XP(NL) 

YY = YP(NL) 

DO 9997 1 = 1, м 

ED = aN e SCALE 

YP(I) = YP(I) * SCALE 

CALL SING(NU,N,XP, YP, XLE, ZLE, TEA, TES, XSING, YSING,CHD) 
CALL TABINT(XP,, YP, XSING, YSING,N) 
NBODY » IE * 1 - IT 

DO 6791 I = 1 , NBODY 

ү|ж1-1 

d - ХР(1) 

Е(1144) = ҮР(1) 

IEP1 = IE + 1 

SLOPT = TES « .75 

DO 438 1 x IEP1 , II 

11 I *1 - IE 

Е(І) ж АӘ(11,1) 

DXI = 1. / 48. 

D = 4. / 3. * (E(I) - .25) 

F(I) = F(IE) + SLOPT * ALOG(D) / D 

L = ІІРІ - I 

E(L) » E(I) 

F(L) = F(IT) * SLOPT * ALOG(D)/D 
WRITE (6,439) 

FORMAT(2X,3H  I,19X, 1HX,19X , 1HY) 
WRITE (6.37) (I;ECI);ECI).I -- 1. 11) 
CALL WRAP(II,JJ,XSING, YSING, E, F,S0,A0, YO) 


C!!II!MAP GRID BACK TO PHYSICAL PLANE AND SHIFT TO QUARTER CHORD 


C 
C 
C 


10 


37 


DO 10 J 22 , JJ 
DO 10 I1 2 1 , II 
Х(1,/-1) ж АӘ(1,1)е»2 - (59(1,1)%Ү(,,1))е%2 


1 - 0.25 


Z(1,J-1) = 2. « АӘ(1,1) « (59(1,1)-Ү9(,,1)) 
JJ = Jd =~ 1 

RETURN 

ҒОКМАТ( 15,2Ғ29.8) 

ЕМО 

SUBROUTINE CLUSTR(DS) 
COMMON/GRID1/X(161,41),2(161,41) 
COMMON/DGRID/DT , IMAX ,KMAX , ITEL, ITEU 
DIMENSION S(41),XP(41), YP(41),R(41) 


THIS SUBROUTINE CLUSTERS A GIVEN X,Z GRID SUCH THAT THE FIRST POINT IS AT 
THE USER-SPECIFIED DISTANCE DNMIN 
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19 


20 


nara — ат 
— ome A) = 


S(K) = SORT((XP(K)-XP(K-1))+*.2+(YP(K)-YP(K-1))0.2) 


1%5(К-1) 


SUMDX = S(KMAX) 

CALL STRTCH(SUMDX,DS ,F1,KMAX, FACTOR) 
WRITE (6,200) I, FACTOR 
R(1) = д. 

DR = DS 

DO 20 K = 2 , KMAX 

R(K) = R(K-1) 4 DR 

DR = DR + FACTOR 

CONTINUE 

RLAST = 1. / R(KMAX) 

DO 3@ K = 2 , KMAX 

R1 » R(K) * RLAST * S(KMAX) 


CIIIII REDISTRIBUTE THE CONSTANT-ETA LINES 


30 
100 


110 
115 


120 
200 


O O O O O 


OO 


CALL TAINT(S,XP,R1,XP1,KMAX, 3, NER, MON) 
Х(1,К) = ХР1 
CALL TAINT(S, YP,R1, YP1,KMAX, 3, NER, MON) 
Z(I,K) = YP1 
CONTINUE 
CONTINUE 
WRITE (6,115) 
DO 110 I = 1 , IMAX 
DX = 212) - Х(1,1) 
ОУ = 2(1,2) - 2(1,1) 
DN » SQRT(DX«DX4DY «DY) 
WRITE(6,120) I , DX , DY , ON 
CONT INUE 
RETURN 
FORMAT (5X , HNORMAL , 1X, 8HDISTANCE,3H AT,4H THE,5H МАСС, / 


1,5H 1,8Х,2НОХ,8Х,2НОҮ,8Х,2НОМ,//) 


FORMAT(15,3F10.5) 

FORMAT(I5,F10.5) 

END 

SUBROUTINE STRTCH(SUMDX,0X1,F1,N1,R) 


THIS SUBROUTINE DETERMINES A GEOMETRIC 

PROGRESSION OF GRID SPACING BETWEEN 1 AND N1 SUH THAT 
SUM8DX) EQUALS SUMOX. THE RATIO BETWEEN SUCCESSIVE 
SPACINGS IS R. 

М = М1 - 1 

К = 1.5 

Е1 = 1.Ё-04 

E2 « 1.Е-04 

00 19 L = 1, 50 

F= (R-1) * 50МОХ - 0Х14(БеэМ-1) 

FP = SUMDX — DXI + FLOAT(N) e R ee (N-1) 

RITER = R - F/ FP 
[F(1.E-02.LT.RITER.AND.RITER.LT. 1.) RITER = 1. 
IF(1..LT.RITER.AND.RITER.LT.100.) RITER=.01 
IF(ABS(R-RITER).LT. R«E!) GO TO 1 

R = RITER 


10 CONTINUE 
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— 


R = 1.0001 

DX1 = DZTOT/FLOAT(N1-1) 

RETURN 

R= RITER 

RETURN 

END 

SUBROUTINE EDDY(DALFA) 
COMMON/FLOW/Q1(161,41),Q2(161,41),Q3(161,41) ,Q4(161,41) 
COMMON /MUTUR/CMU( 161,41) 

COMMON/SK INCF/CF (161) 
COMMON/DGRID/DT , IMAX, KMAX , ITEL, ITEU 
COMMON/PAR/GAMMA , REYREF , ALFA, ALFA1, REDFRE , AMINF , ALFAT 
COMMON/GR ID1/X(161,41),2(161,41) 

DIMENSION TIN(41),TOUT(41), Y(41) 


INITIALIZE VISCOSITY EVERYWHERE 
FACT! = OT + AMINF / REYREF 
СМОМАХ = 100. » FACT! / DT 
DO 1 K = 1 , KMAX 
DO 1 I= 1 , IMAX 
CMU(I,K) = FACTI 
THIS SUBROUTINE COMPUTES THE EDOY VISCOSITY BASED ON THE 
BALOWIN-LOMAX TWO LAYER MODEL 


1 = 2 , ІМАХ - 1 
UDIF = 0. 

0.1Е-06 
. 1E-06 


101(1:1) 


COMPUTE TAU AT THE WALL 


МЕТ « 1.4(02(1,2)/01(1,2) - аа аиа 
ҮЕТ » 1.4(03(1,2)/01(1,2) - 03(1,1)/01(1.1) 
XXI = 21 - ы 

ZXI * Z(I*1,1) - Z(I-1,1) 

ХЕТ = 4. , Х(1,2) - 3. « X(I,1) - X(1,3) 
СЕТ = 4. • 2(1,2) - 3. * Z(I,1) - Z(I,3) 
ХХІ = .5 ® ХХ] 

ZXI = .5 » ZXI 

ХЕТ = .5 • ХЕТ 

LEN a, es ЛЕТ 

YAC = 1. / (XXI * ZET Е * XET) 

OMEGA = (UET * XXI – VET = ZXI ) * YAC 
TWALL = AMINF * OMEGA / REYREF 

CF(I) = 2. * TWALL / (AMINF *e2) 


FACT = SQRT(Q1(1,1) * ABS(TWALL))*REYREF/(26.*AMINF) 
DO 10 K = 2 , KMAX-1 





UX] = (Q2(1+1,K)/Q1(1+1,K) - Q2(I-1,K)/Q1(1-1,K)) 
VXI = (Q3(1+1,K)/Q1(1+1,K)-Q3(I-1,K)/Q1(I-1,K 
UET = (Q2(1,K+1)/Q1 foede 701(1,К-1 
VET » (03(1,К-1)/01(1,К-11-03(1,К-1)/01(1,К-1 
XXI sw Х(141,К) - Х(1-1,К 
ZXI = Z(I+1,K) — Z EE, 
XET = X(I,K4*1) - Х(1,К-1 
СЕТ = 7(1,К+1) - 2(1,К-1) 
ТАС = 
А 


TA EE е 7ЕТ - ZXI * XET) 


e ABS(UETeXXI*VET*eZXI-UXI*XET-VXIeZET) + YAC 
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31 


ОООО 


10 


20 


UDIF = SQRT(Q2(I,K)**24Q3(1,K)**2)/Q1(I,K) - UWALL 
IF(ABS(UDIF).GT.UDIFMAX) UDIFMAX = ABS(UDIF) 

Ү(К) « 5О8Т((Х(1.К)-Х(1,К-1))6424(2(1,К7-2(1,К-1))в2)4Ү(К-1) 
Ғ ж Ү(К) * OMEGA 

e GO TO 31 

IF(I.GT.ITEL. ANO. I. LT. ITEU) F = F * (1. — EXP(—Y(K)eFACT)) 
CONT INUE 


MODIFIED TURBULENCE MODEL APPLY FOR SPECIFIED RANGE OF ANGLES WHERE 
FY IS USED TO FIND THE SECOND PEAK VALUE OF F FUNCTION 


IF(ALFA. LT.ALFAI.AND.DALFA.GE.0.) THEN 
FY = Е » Y(K) 
IF(FY.GT.FYMAX) THEN 
FYMAX FY 
FMAX 
YMAX 
ENDIF 
ENDIF 
IF(ALFA.GE.ALFAI.OR.DALFA.LT.@.) THEN 
IF(F.GT.FMAX) THEN 
FMAX = F 
YMAX = Y(K) 
ENDIF 
ENDIF 
FCT = Y(K) * FACT 
IF(FCT.GT.20.) FCT = 20. 
FCT = ABS(FCT) 
EL = .4 « Y(K) » (1. - EXP(-FCT)) 
TIN(K) 2 Q1(I,K) * EL * EL * OMEGA 
TIN(K) = ABS(TIN(K)) 
CONT INUE 
KSWTCH ж д 
FWAKE = YMAX » ҒМАХ 
FI ж 0.25 « ҮМАХ «e UDIF #2 / FMAX 
IF(F1.LT.FWAKE) FWAKE = F1 
DO 20K 22 , KMAX - 1 
FKLEB = 0. 
IF(ABS(Y(K)/YMAX).LT.1.E+04) THEN 
FKLEB = 1. / (1. + 5.5 « (0.3 » Y(K)/YMAX) e. 6) 
END IF 
TOUT(K) = .0168 * 1.6 » 01(1,К) » РИАКЕ « РКЕЕВ 
TOUT(K) = ABS(TOUT(K)) 
E 0) СО ТО 20 
IF(TIN(K).GT.TOUT(K)) KSWTCH = K – 1 
CONT INUE 


Е 
Y(K) 


CITIIITOTAL VISCOSITY IS SUM OF LAMINAR AND INNER/OUTER LAYER AS APPROPRIATE 


3e 
100 


DO 30 К = 2 , KMAX - 1 

IF(K.LE.KSWTCH) THEN 

Ee) = ОТ ® (AMINF/REYREF * ABS(TIN(K))) 
LSE 

CMU(I,K) = DT © (AMINF / REYREF + ABS(TOUT(K))) 

END IF 

CONT INUE 

CONTINUE 

RETURN 

END 

SUBROUT INE RESI(OMEGA) 

COMMON/PERTR/DQ1 (161,41),DQ2(161,41),0Q3(161,41),004(161,41) 

COMMON/GR ID1/X(161,41),2(161,41) 
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COMMON/DGR ID/DT , IMAX, KMAX , ITEL, ITEU 
COMMON/FLOW/Q1(161,41),02(161,41),Q3(161,41),04(161,41) 
COMMON/PAR/GAMMA , REYREF , ALFA, ALFAT, REDFRE, AMINF , ALFAI 
DIMENSION RHS(161,4) 

e = OMEGA « Z(1,K) 

YTAU( I,K) = — OMEGA « X(I,K) 

THIS SUBROUTINE COMPUTES THE RESIDUAL ON THE RIGHT HAND 

SIDE ARISING FROM THE EULER- PART OF THE GOVERNING EQUATIONS 


FLUX TERMS WITHIN THE XI- DERIVATIVE 

DO 100 K = 2 , KMAX — 1 

DO 10 1 = 1 , IMAX 

ОСОМ = (Q2(1,K)/01(1,K) ) ® (2(1,к+1)-2(1,к-1)) 
1 — (Q3(1,K)/Q1(I,K)) © (X(I,K*1)-X(I,K-1)) 

ОСОМ жз 0.25 e DT * UCON 

XIT = — XTAU(I,K) *(Z(I,K*1)-Z(I,K-1)) 

1 + УТАЦ(Т,К) « (Х(1,К41) - Х(1,К-1)) 

ХІТ жз ХІТ * DT * 0.25 

UCON = UCON + XIT 

RHS(I,1) s UCON e Q1(1,K) 

R = 1. / QI(I,K) 

P = (GAWA-1.) * (Q4(I,K) - .5 * R *(Q2(1,K)0.2+ 
1 Q3(1,K)**2)) 


ВН5(1,2) = 02(1,K) + UCON + P + DT «+. 9.25 e dee. - Z(I,K-1)) 
RHS(1,3) = Q3(1,K) « (ОМ -Р « 0Т «09.25 « (Х(1,К-1)-Х(1,К-1)) 
RHS(1,4) = UCON * (Q4(I,K)+P) – ХІТ • Р 


10 CONTINUE 
DO 11 1*€ 2 , IMAX — 1 
001(1,К) = DQ1(I,K) — RHS(I+1,1) + RHS 
DQ2 nr = DQ2(I,K) - RHS(I+1,2) + RHS 
003(1,К) = DQ3(1,K) - RHS(I+1,3) + RHS 
11 DQ4(1,K) « 004(1,К) - RHS(I+1,4) + RHS( 
100 CONTINUE 


1-1, 1 
1-1,2 
1-1,3 
1-1,4) 


FLUX TERMS WITHIN THE ETA- DERIVATIVE 


DO 200 1 2 , IMAX — 1 
DO 20 K = 1 , KMAX 
VCON = Иа ) а ER 
4(03(1,К)/01(1,К . (Х(141,К)-Х(1-1,К 

VCON в УСОМ в 0.25 » DT 

ETAT = -XTAU(I,K) * (Z(I-1,K)-Z(I*1,K)) - YTAU(I,K)e 
1 (Х(141,К)-Х(1-1,К)) 
ETAT ж ЕТАТ » 0.25 « ОТ 

УСОМ = УСОМ + ЕТАТ 

RHS(K,1) = VCON e RT 

Р ж (САММА-1.) , (04(1,К) - 9.5 s«(Q2(I,K)**24Q3(I,K)**2)/Q1(I, K)) 
к = VCON e Sd +P * ОТ * .25 e (Z(I-1,K)-Z(I*1,K)) 


1 


RHS(K,3) = VCON * Q3(I,K) + P = DT e .25 * (Х(1+1,К) — X(I-1,K)) 

RHS(K,4) = VCON « (Q4(I,K)*P) — ETAT « P 
20 CONTINUE 

DO 21 К = 2 , KMAX - 1 

001(1,К) » OQI(I,K) — RHS(K+1,1) + ЕН5(К-1,1 

002 n = 002 LK) - Е * RHS(K-1,2 

003(1,К) = DQ3(1,K) — RHS(K+1,3) + RHS(K-1,3 
21 DQ4(I,K) = DQ4(1,K) – RHS(K+1,4) + RHS(K-1,4 
200 CONTINUE 

RETURN 

ENO 

SUBROUTINE ROTGRID(X,Z, IMAX,KMAX,DALFA) 
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С ROTATE GRID IN THE CLOCKWISE DIRECTION BY AN AMOUNT DALFA 

DIMENSION X(161,41),2(161,41) 
CA = COS(DALFA) 
SA « - SIN(DALFA) 
DO 20 K = 1 , KMAX 
DO 20 I = 1 , IMAX 
X1 = X(I,K) 
21 = Z(I,K) 
2006) = X1 * СА — Z1 * SA 

20 Z(I,K) = Z1 * CA + X1 e SA 
RETURN 
END 
SUBROUTINE CPPLOT(I1,12,FMACH,X,Y,CP) 


THIS SUBROUTINE PLOTS CP AT EQUAL INTERVALS IN THE MAPPED PLANE 


O O O 


DIMENSION KODE(4),LINE(90),X(161),Y(161),CP(161) 
DATA KODE/1H ,1H+,1HI,1He/ 
WRITE ( 6 , 2) 
2 FORMAT(5@H@PLOT OF CP AT EQUAL INTERVALS IN THE MAPPED PLANE/ 
1 10H90 X , 19H X/C ,10H CPL ,10H CPU ) 
CPO = (1. + .2 e FMACH +22) e. 3.5 -1. 
CP@ - CPO / ( .7 * FMACH #22) 
KO = 30. « CPO + 4.5 
ІМІМ ж (12-11)/2 + 11 
ILOW = 2 « IMIN 
CHO=X( 11) - X( IMIN) 
DO 12 I 1 , 90 
12 LINE(I) KODE(1) 
DO 34 I IMIN , 12 
К ж 30. (СРд - СР(1)) + 4.5 
K1 = 30. * (CP@ -CP(ILOW-I)) + 4.5 
IF(K.LT.1) K = 1 
т K1 = 1 
IF(K.GT.90) K » 90 
IF(K1 .GT. 98) K1 = 90 
(ІМЕ(кӨ) = KODE(3) 
АЕО = KODE(2) 
LINE(K1) = E 
XOC = (X(I) - X(IMIN)) / CHD 
WRITE (6,610) Шо осе уын 
Re EE ж КООЕ(1) 
34 LINE(K) = KODE(1) 
RETURN 
610 FORMAT(4F10.4,90A1) 
END 


en il 


/toF 
NACA 0012 AIRFOIL, RUN: 
159 41 ‚0925 ә. 15.000 19.09 0.90 0.151 "5999 
МО. ОР 5ТЕР5 
4509. 
REYNOLDS NUMBER IN MILLIONS, DISTANCE OF FIRST POINT OFF THE WALL 
‚345 . 00005 
TSTART 
9.0 
FULOUT 
-1.9 
RESTART ,PITCH,PLUNGE ,OUTPUT OPTIONS SPECIFIED IN THE NEXT CARD 
TRUE TRUE FALSE 
FNU FNL FSYM 


lll 


23. 


:01155 
.01894 
.02615 
2039955 
.04200 
.04683 
.05345 
.05737 
. 05941 
. 05962 
. 05978 
. 05992 
. 05999 
. 06000 
. 05999 
. 05992 
. 05980 
.05965 
. 05947 
. 05803 
. 05581 
. 95294 
. 04952 
. 04563 
. 04137 
‚03664 
. 03161 
‚02623 
‚02055 
‚01448 
. 00807 
. 00126 


ІШ? 


APPENDIX B 
NOTES ON USE OF THE NAVIER-STOKES SOLVER 


l. |. JOB CARDS 

The JCL options are selected bv removing the “comment” designator ("*.”) at the 
beginning of the applicable lines. The options available are: 

өг save tie current solution, Values of TIME, Q1, Q2, Q3, and Q4 are saved 
(logical unit 08) for subsequent restart. Activate the two lines referencing 
ЧЕ. 

e Create pressure coefficient data file. Data is accumulated through the runs to 
be accessed and read by a separate program. Activate the “OLDCP” statements 
to access and add to previously stored data. Activate " NEWCP" statements to 
store current cumulative data. A new version 1s created each time, so the files 
must be purged periodically. 

e Start from a stored solution. The above values are read (logical unit 07) and 
Iteration continued from that point. Activate the two lines referencing 
SOW SIN’. 

e Creat PLOT3D files. Conversion from Cray to VAX binary is handled and 
properly formatted “Q” and “XYZ” files are created for use with the PLOT3D 
graphics program. Activate the Q and XYZ “ASSIGN” statements and the 
“SENDVAX” and “ACCESS” statements before the “EXIT” card. 

өе Create “troubleshooting” PLOT3D files. If the solution “blows up” and the 
appropriate “WRITE” statements are not commented out in the main program, 
Q and XYZ files are created to investigate the status of the solution at the last 
“successful” iteration. Activate the “ASSIGN” statements and the “ACCESS” 
and “SENDVAX” statements after the “EXIT” card. 

e Use job chaining. The FETCH, REWIND, and SUBMIT statements are used 
to call up another program when current run is completed. 

In addition, if it is desired to save more than one solution, the PDN (5 digits) 
and ID may be changed. Old data sets must be purged from Cray. This is 
accomplished by placing an “AUDIT.” card after the account card when running a 
program or by running “AUDIT.JCL” to obtain a listing of current data sets. Then, 
В, ма SRILTICL to create a JCL to delete Cray PDN s. It is then only 
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necessary to precede this JCL by a JOB card and ACCOUNT card and Tune 


program as usual. Consult the ACF Cray Users’ Guide for detailed information on job 


control cards. 


p 


MAIN PROGRAM 


Certain changes may be made within the main program. These changes affect 


the program execution or change the output. 


3. 


The frequency of steady-state output may be changed bv varying the value in 
the first statement after the comment “FOR STEADY STATE OU TPU iii 
TRE РОТЕ. 

The interval for exiting the program (in order to save a solution and generate 
PLOT3D files) may be changed under *MULTIPLY ХхСРОС Т Ун 
example: NEXIT = 24 * NCPOUT means the program will exit automatically 
four times a cycle, or every 24 times the normal printed output is generated. 
Velocity profile information may be output if desired (as when a permanent 
record of a converged solution is desired) bv activating the WRITE statement 
just before the CONTINUE statement numbered 4000. This outputs the 
indices of the X and Z coordinates of each grid point with the corresponding 
values of pu, pv, eddy viscosity, total velocity (/u* + v*), and distance normal 
to the wall. 

At the beginning of subroutine SLPS, the value of the implicit damping 
coefficient may be adjusted by changing the multiple of WW. 

At the end of subroutine SLPS, the frequency with which the residuals are 
output may be set by changing the value under “SELECT ІНЕ ТУТЕКУА 
At least every ten steps 1s recommended. 

Other output values may be turned on and off as well. “UNWRAPPED 
COORDINATES IN THE TRANSFORMED PLANE” is in subroutine 
“WRAP”. Airfoil coordinates are in “AIRFOL”. Grid spacing and “NORMAL 
DISTANCE AT THE WALL arein Cleo 


DATA CARDS 


Most of the “tuning” of the program is done with the data cards: 


Ist Line - Name of airfoil and other identifying information. Eighty characters 


available. 
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2nd Line - (1) and (2) The first two values, IMAX and KMAX (format 15), аге 
the number of X and Z coordinate locations to be used. These remain at 159 
and 41 for the 161 X41 C-grid used at present. The remainder of the line 
contains seven values in format F10.0. (3) DT: size of the time step. This is 
automatically set to 1.0 within the program when the reduced frequency is less 
than or equal to 0.001. In this case the program uses the local time-step option 
(relaxation). For dynamic stall DT = 0.005 is recommended, at angles below 
20 degrees, smaller at higher angles. (4) WW: explicit artificial viscosity term. 
Normally 2-5, with about 2-3 recommended for static cases, 5 for dynamic stall. 
Higher numbers have greater effect on solution. (5) ALFA: mean angle of 
attack. (6) ALFAI: amplitude of oscillation. (7) ALFAI: angle below which 
a modified turbulence model is used (upstroke only) to compute eddy viscosity 
(Baldwin-Lomax model). Normally set to 19 degrees for dynaniuc stall. May 
ШРСИФТАрШаР СІР olution. (Ss) REDERE: reduced frequency. (9) AMINF: 
Mach number. 

3rd Line - “NO. OF STEPS”--not read. 

4th Line - FNSTP: number of time steps to be done on the present run (format 
F 10.4). 

3th Line - “REYNOLDS NUMBER...”--not read. 

6th Line - (1) REYREF: the Revnolds number in nullions (format F10.4). A 
negative value means inviscid flow. (2) DNMIN: distance of the first point off 
the wall (format F10.4). For Reynolds numbers up to 3 nullion 0.00005 should 
be used normally. Stabilitv of the solution may be improved bv increasing this 
value in some cases (1e., high AOA steady angle of attack). 

7th Line - TSTART --not read. 

Sth Line - TSTART (format F10.4): time calculations have been advanced up 
to the previous run. When a negative value is used, TSTART is read from 
logical unit 08 (see JCL comments). Then normally use 0.0 for initial runs and 
—1.0 for restarts. Must use 0.0 for first dvnamic run from converged steadv- 
state solution. 

9th Line - "FULOUT "--not read. 

lOth Line - FULOUT (format F10.4): —1.0 means no plotting files will be 
generated. Set 0.0 to begin full output, then set 1.0 to continue. When using 
full output, the appropriate job cards must be activated. 

Lith Line - RESTARTI, PITCH- -not read. 
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e 12th Line - RSTRT, PITCH, PLUNGE (format E5) ІГ ЕЗТЕТ 15 50 104 
values of TIME, Ql, Q2. Q3, and Q4 are read to continue iteration ГГГС A 
set for dynamic stall and indicates change in angle of attack. PLUNGE will not 
be set for present purposes. It indicates up and down motion of the airfoil. 


e 13th Line - The remaining lines define airfoil geometry and for the present are 
set for the NACAOO12 airfoil. 


4, ADDITIONAL NOTES 

For high angle of attack (separated flow), or if otherwise desired for time 
accurate steady-state calculations, use reduced frequency of 0.002 and ALFA] 
(amplitude) = 0. Set PITCH = FALSE and DT = 0.005, as for dynamic stall. For 
stability WW should be set to 5 (or even higher, up to 10) and ALFAI at or below the 
minimum angle to use the original Baldwin-Lomax turbulence model. The distance of 
the first point off the wall mav also be increased ( ^ 0.0001), which has the effect of 
coarsening the grid for the troublesome fine-mesh leading-edge area. 

When doing dvnamic stall simulations where the AOA goes to the 20-25 degrees 
range, use the DT = 0.005 at lower angles for computation time, but reduce this value 
when restarting going into the high angle portion of the cycle. 

The values of DRMAX, DUMA(X, etc. may be useful when a solution "blows 
up”. The IR and KR values give the indices of the X - Z location on the grid where 
density changed the fastest (IR — $0 is the leading edge). Normally. this will be near 
the leading edge. DT should then be reduced. 

Revnolds numbers of 109 to 10 x 10 should show no effect on sensitivity of 
solution for distance of the first point off the wall of .00005, since this places several 
points in the boundary laver. 

Mach numbers of .1-.5 should not alter calculation time significantly, but 


perhaps twice as many iterations may be required for convergence in the .5-.8 range. 
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